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I.  INTRODUCTION 

The  original  Douglas  Neumann  Cascade  Program,  documented  in  a  report 
by  Joseph  P.  Geising  [1]*,  solves  an  approximate  potential  solution 
for  the  flow  through  two-dimensional  linear  cascades.  This  approximate 
solution  approaches  the  exact  solution  as  the  number  of  points  describing 
the  body  surface  approaches  infinity.  The  original  program  can  handle 
one  or  more  cascades,  the  bodies  can  be  of  any  arbitrary  geometry,  and 
the  bodies  can  be  lifting  or  nonlifting. 

The  Applied  Research  Laboratory  of  The  Pennsylvania  State  University 
modified  the  Douglas  Neumann  program  so  that  flow  inward  through  a  single 
radial  flow  or  circular  cascade  could  be  analyzed.  This  modified  program  was 
developed  under  sponsorship  by  the  U.S.  Bureau  of  Reclamation  to  enable  the 
flow  through  wicket  gates  of  a  hydraulic  turbine  to  be  analyzed  [2],  A  con¬ 
formal  transformation  was  incorporated  into  the  program  which  transformed  the 
circular  cascade  into  a  linear  cascade.  The  transformed  cascade  could  then  be 
analyzed  in  the  conventional  manner  and  the  solution  then  transformed 
back  to  the  real  plane.  In  Reference  [2],  average  downstream  flow  angles 
calculated  with  the  modified  program  are  compared  to  experimental  values 
obtained  in  an  air  test  facility.  The  agreement  between  the  experimental 
and  calculated  flow  angles  is  found  to  be  quite  good. 

When  the  work  of  [2]  was  performed,  much  of  the  flexibility  of  the 
original  Douglas  Neumann  Cascade  program  was  lost,  and  also,  both  linear 
and  circular  cascades  could  not  be  analyzed  with  a  single  program.  This 
report  documents  a  rewritten  version  of  the  Douglas  Neumann  program  which 
can  handle  both  linear  and  circular  cascades,  multiple  cascades  of  both 
types  can  be  analyzed,  and  the  input  parameters  which  define  the  flow 


*Numbers  in  brackets  refer  to  references  at  the  end  of  the  report. 
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conditions  through  the  cascades  can  be  of  several  different  types  as  in 
the  original  program.  In  addition,  the  new  program  can  analyze  circular 
cascades  with  the  radial  flow  direction  inward  or  outward.  Thus,  for 
example,  the  flow  through  stay  vanes  and  wicket  gates  of  a  pump-turbine 
can  be  analyzed  in  both  the  pump  and  turbine  mode. 

With  the  many  different  types  of  problems  to  be  handled  by  a 
single  program,  it  was  thought  that  the  program  would  be  more  orderly 
if  it  was  rewritten  following  the  original  technique,  rather  than 
inserting  a  lot  of  new  sections  into  the  old  program.  The  program  is 
now  structured  so  that  data  are  read  in,  manipulated,  calculations 
are  performed  etc.  in  progressive  steps  through  the  program.  There  are 
numerous  sections  of  the  program  which  will  not  apply  to  a  particular 
problem,  thus  the  program  skips  these  sections  and  continues  on  with 
the  next  section  which  applies.  The  program  does  not  double  back  over 
itself,  except  in  small  iteration  loops,  and,  thus,  the  logic  should  be 
easy  to  follow.  The  functions  of  the  various  sections  of  the  new  program 
are  also  documented  with  comment  statements  throughout  the  program. 

In  the  Douglas  Neumann  Cascade  program,  the  flow  can  be  analyzed 
at  points  off  the  body  in  addition  to  obtaining  the  flow  on  the  body 
surface  and  the  overall  cascade  performance.  In  many  circumstances  it 
may  be  desirable  to  analyze  a  large  number  of  off-body  points  which  are 
systematically  spaced  through  the  cascades.  The  rewritten  version  of  the 
program  will  generate  off-body  point  coordinates  given  an  initial  set  of 
coordinates  and  the  desired  spacing.  In  the  new  version,  the  off-body 
point  coordinates  are  also  stored  in  arrays  separate  from  the  body 
coordinates  so  that  the  number  of  off-body  points  need  not  be  considered 
when  selecting  the  number  of  points  to  represent  the  airfoils.  Two 
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large  arrays  have  been  specified  for  storing  off-body  point  coordinates 
and  the  separate  storage  makes  changing  the  dimensions  of  the  arrays  an 
easy  task,  if  it  becomes  necessary. 

The  use  of  auxiliary  and  core  storage  has  also  been  greatly  modified 
in  the  new  version  of  the  program.  The  original  Douglas  Neumann  program 
used  tapes  to  store  data  and  transfer  the  data  from  subroutine  to  subroutine. 
This  was  probably  done  so  that  only  one  subroutine  needed  to  be  in  the 
core  at  one  time,  and  thus  the  program  could  be  run  on  smaller  computers. 

With  today's  larger  computers,  COMMON  storage  can  be  used  for  transferring 
data  between  subroutines,  and  the  repetitious  reading  of  data  from  auxiliary 
storage  eliminated.  Auxiliary  storage  is  currently  only  used  to  store 
the  original  body  coordinates  and  the  very  large  matrices  representing 
the  system  of  equations  which  are  solved  simultaneously. 

For  very  large  problems,  the  matrix  representing  the  system  of 
equations  to  be  solved  is  too  large  to  fit  in  the  core  of  most  computers. 

The  original  Douglas  Neumann  program  utilized  an  iterative  technique  and 
auxiliary  storage  to  solve  the  equations  for  large  problems.  The  iterative 
approach  required  reading  the  equations  from  a  file  for  each  iteration.  Since 
reading  of  disc  files  is  time  consuming,  a  new  procedure  for  solving  the 
equations  was  developed  during  the  current  rewriting  of  the  program.  Both 
the  new  procedure  and  the  iterative  procedure  are  available  with  the  rewritten 
version  of  the  program.  The  new  procedure  was  shown  in  a  sample  run  to  signifi¬ 
cantly  reduce  the  computation  time. 

Thus  far,  this  introduction  has  given  a  general  description  of  the 
background  and  development  of  the  cascade  program.  The  major  changes  made 
in  the  rewritten  version  of  the  Douglas  Neumann  cascade  program  have  also 
been  briefly  described.  In  the  remainder  of  this  report,  the  theory  used 


i 
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in  the  program,  a  description  of  the  program  and  the  input  and  output 
parameters,  and  several  check  cases  will  be  presented.  A  users  guide 
giving  detailed  instructions  for  running  the  program  and  a  sample  run 
can  be  found  in  Appendix  A.  Appendix  B  discusses  the  data  storage  and 
solution  procedures,  and  Appendix  C  describes  the  coordinate  data  input 
order  and  the  necessary  reordering  performed  in  the  program. 


II.  THEORY 
2 . 1  Basic  Equations 

The  basic  concept  of  the  Douglas  Neumann  cascade  program  is  to  apply 
a  distribution  of  sources  on  the  surface  of  the  bodies  such  that  the 
combination  of  the  source  distribution  and  the  onset  flow  satisfy  the 
boundary  condition  on  the  blade  surface.  For  the  usual  condition  of  no 
flow  normal  to  the  surface,  the  following  equation  should  be  satisfied 
everywhere  on  the  blade  surface: 

~Vm  •  Ti(s)  =  |  a(q)  A(q , s)  dq  (1) 

body 


where 


V  =  the  onset  flow 
00 

~n(s)  =  unit  vector  normal  to  the  surface  at  point  s 

o(q)  =  source  strength  at  point  q 

-A  v 


A(q,s)  =  V(q,s)  •  n(s) 

V(q , s)  =  velocity  at  point  s  due  to  an  infinite  array  of  unit  sources 
at  point  q  on  the  bodies. 

An  expression  for  V(q,s)  can  be  obtained  by  starting  with  the  equation 
for  the  velocity  at  a  point  due  to  a  single  source  and  writing  V(q,s) 
as  an  infinite  sum  of  the  velocities  induced  by  an  infinite  array  of 
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sources.  The  infinite  sum  can  then  be  eliminated  by  manipulating  the 
equation  and  introducing  the  hyperbolic  functions.  The  following  equation 
is  obtained  for  V(q,s): 

~  t*(a)  "  C(q)3  ’  (2) 

where 

z(s)  ■  complex  coordinate  at  point  s  =  x  +  iy 
4(q)  =  complex  coordinate  at  point  q  =  £  +  irj 
The  expression  for  the  velocity  induced  by  an  infinite  array  of  vortices 
is  the  same  as  Equation  (2),  but  with  the  velocity  vector  rotated  90°. 


V(q,s)  =  V  -  iV  = 

^  source  x  y 


V(q, s) 


vortex 


~  tz(s)  -  C(q)] 


The  technique  employed  in  the  Douglas  Neumann  program  to  solve 
Equation  (1)  is  to;  (1)  break  the  surface  of  the  bodies  into  small  segments 
or  elements,  (2)  approximate  the  source  distribution  as  a  constant  on  each 
element,  and  (3)  satisfy  the  boundary  condition  only  at  the  element  midpoint. 
With  these  approximations.  Equation  (1)  can  be  written  as: 

J*  »  f 

vooi  ‘  ni  =  l  ° k  Mq)  d(q)  (4> 

3  3  k=i  K  ,  >  ,  3 

element  k 

As  seen  in  Equation  (4),  the  independent  variable  s  has  been  eliminated, 
because  Equation  (4)  is  evaluated  only  at  the  midpoints  of  the  elements. 

Thus,  A(q,s)  becomes  Aj (q)  and  can  be  written  as: 

Aj(q)  =  Vj(q)  •  n.  =  coth  {gj  [Zj  -  4(q)  ]  }j  •  *n.  (5) 

where 

z.  =  the  complex  coordinate  of  the  midpoint  of  element  j. 
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The  integral  in  Equation  (4)  can  be  evaluated  by  expressing  the  surface 
coordinate  q  in  terms  of  the  global  coordinates  of  the  body  in  the  cascade. 
With  the  aid  of  Figure  1  we  find  that: 


£(<j)  =  (Cc  +  i  nQ)  +  q (cos  £>k  +  i  sin  g^)  =  Cq  +  qei,Jk  (6) 

k 

Equation  (6)  can  then  be  differentiated  to  yield  an  expression  for  dq, 
which  is  needed  to  perforin  the  integration. 


dq  =  d£e“iek  (7) 

Substituting  Equations  (5)  and  (7)  into  the  integral  and  expressing  the 
integration  limits  as  the  end  points  of  element  k  as  defined  in  Figure  1 
yields: 

f  2k 

J  Aj(q)  dq  =  "ru  •  j  •"  coth  (Zj  -  O]  d£  e-1^  .  (8) 

element  k  r 

1k 


Carrying  out  the  integration  in  Equation  (8),  provides  an  expression  for 
A  ,  the  normal  velocity  induced  at  the  midpoint  of  element  j  due  to  the 
source  distribution  on  element  k  of  constant  unit  strength. 


jk 


-iSk 

Aj  (q)  dq  =  —  In 


element  k 


sinh  ^  (z.  -  Clk) 
sinh  gj  (z.  -  C2k) 


(9) 


Substituting  A  as  defined  in  Equation  (9)  into  Equation  (4)  yields  the 
j  k 

basic  equation  used  in  the  Douglas  Neumann  program. 

X  V  ■  %  •"”)  (10) 

Equation  (10)  is  the  approximate  form  of  Equation  (1)  and  it  must  be 
satisfied  for  each  element  of  the  cascades.  Equation  (10)  can  be  written 
for  each  element  yielding  the  following  matrix  equation: 


-16- 


26  June  1981 
AMY : cag 


(11) 


The  set  of  equations  represented  by  Equation  (11)  are  solved  simultaneously 
in  the  Douglas  Neumann  program  for  the  unknown  source  strength  of  each  element. 
As  seen  by  examining  Equation  (9),  the  terms  of  the  A  matrix  are  only  functions 
of  element  coordinates  and,  thus,  can  be  calculated  from  the  cascade  geometries 
which  are  input  to  the  program.  The  onset  flows,  appearing  in  the  matrix  on 
the  right  hand  side  of  Equation  (11),  will  be  discussed  in  the  following 
section. 

2. 2  Basic  Solutions  and  Combination  Equations 

Equation  (11)  can  be  solved  for  any  given  onset  flow  V^.  However, 
for  airfoils  in  a  cascade,  an  addition  boundary  condition,  must  also  be 
satisfied.  Circulation,  which  arises  around  the  airfoils,  enables  the 
Kutta  condition  to  be  met.  For  cascades,  this  circulation  induces  an 
upwash  and  downwash  velocity  far  upstream  and  downstream,  respectively, 
relative  to  the  mean  cascade  velocity.  These  induced  velocities  are 
illustrated  in  Figure  2.  Both  the  circulation  and  the  mean  velocity  must 
be  accounted  for  in  Equation  (11).  However,  for  the  typical  cascade 
problem,  the  magnitude  of  the  circulation  is  usually  unknown  and  the  mean 
flow  angle  may  not  be  given.  The  operating  parameter  which  may  be  given 
for  a  cascade  analysis  will  be  either: 

1.  cij,  the  inlet  flow  angle 

2.  Aa,  the  angle  the  flow  is  turned  through  the  cascade  =  ot  -  n_ 

I  E 
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3.  a,  the  mean  flow  angle 

or  4.  CL>  the  total  lift  coefficient  of  all  the  cascades. 

The  angles  aT,  a  and  a_  are  shown  in  Figure  2. 

1  L 

In  the  Douglas  Neumann  program,  several  basic  solutions  are  calculated, 
each  of  which  satisfy  the  condition  of  zero  velocity  normal  to  the  body 
surfaces.  Since  the  governing  equation  is  linear,  any  number  of  solutions 
can  be  scaled  and  added  without  violating  the  governing  equation  or  the 
surface  boundary  condition.  The  technique  employed  by  the  Douglas  Neumann 
program  is  to  combine  the  basic  solutions  such  that  the  cascade  operating 
parameter  and  Kutta  condition  are  satisfied. 

The  basic  solutions  which  are  used  are  calculated  for  each  of  the 
following  onset  flows: 

1.  A  uniform  flow  normal  to  the  cascade(s),  i.e.,  a  =  0 

2.  A  uniform  flow  parallel  to  the  cascade(s),  i.e.,  a  =  90° 

3.  A  nonuniform  flow  created  by  the  circulation  on  the 
airfoils  of  cascade  1. 

4.  A  nonuniform  flow,  like  flow  3,  for  each  additional  cascade 
with  circulation. 

The  terms  -V  .  •  n.  in  Equation  (11)  for  onset  flows  1  and  2  are 
°°J  3 

simply  sin(Bj)  and  -cos(Bj),  where  8  is  the  angle  of  the  element  as 
previously  illustrated  in  Figure  1.  For  the  nonuniform  onset  flows, 
the  effects  of  all  the  elements  of  an  airfoil  in  one  cascade  on  a  particular 
element  must  be  summed  to  obtain  :he  velocity  induced  by  the  circulation. 
Since  the  only  difference  in  the  velocity  induced  by  a  vortex  or  a  source 
distribution  is  that  the  velocity  vector  is  rotated  90  degrees,  the  calcula¬ 
tion  performed  in  the  program  to  evaluate  Equation  (9)  can  also  be  used  to 
calculate  the  onset  flows  created  by  the  circulation  of  the  different 
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cascades.  The  normal  component  of  the  onset  flow  for  element  j,  due 
to  circulation  on  the  airfoils  in  cascade  m  can  be  written  as: 


.  r  n 

^  =  —  y2  A* 

Vco-i  nj  g  Z. 


m  k=n. 


where  T  =  the  circulation  on  cascade  n 
m 

S  =  the  circumference  of  an  airfoil  in  cascade  m 
m 

=  the  normal  component  of  the  velocity  induced  at  the  midpoint 
of  element  j  by  a  distribution  of  constant  unit  strength 
vortices  on  element  k.  (Since  the  dot  produce  in  Equation  (9) 
is  performed  by  complex  number  multiplication,  A!,  is  obtained 

JK 

from  the  real  part  of  the  resulting  number  while  A  ^  is  the 
imaginary  part.) 

n^  and  n^  =  the  indices  of  the  first  and  last  elements  of  the  airfoils 
in  the  m  cascade. 

The  value  of  in  Equation  (12)  is  unknown  at  the  time  the  basic 
solutions  are  calculated.  Thus,  the  basic  solutions  for  the  onset  flows 
due  to  circulation  are  obtained  with  T  =1.  The  real  magnitude  of  T 
is  determined  when  the  basic  solutions  are  combined  to  satisfy  the  Kutta 
condition. 

Under  most  circumstances,  the  Kutta  condition  is  violated  by  each 
of  the  basic  solutions.  An  approximate  measure  of  the  amount  each  solution 
violates  the  Kutta  condition  can  be  obtained  by  computing  the  difference  in 
the  tangential  velocity  above  and  below  the  trailing  edge  (i.e.  at  the 
midpoint  of  the  first  and  last  element  of  each  cascade).  The  basic 
solutions  must  be  combined  such  that  the  cascade  operating  parameter 


is  met  and  the  combined  violations  of  the  Kutta  condition  are  zero.  The 
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following  equations  express  these  conditions  and  provide  a  means  for 

evaluating  the  unknown  values  of  T  . 

m 

AV01  U  cos  a  +  U  sin  a  +  AV.^1^  +  AV^I^  +  . . .  =  0 

AV02  U  cos  a  +  AV902  U  sin  a  +  AV2iri  +  AV22r2  +  ...  =  0 

•  •  •  • 

(13) 

•  «  •  • 

•  •  •  • 

AV0NC  u  C0S  a  +  AV90NC  u  sin  a  +  AVNC,iri  +  AVNC,2r2  +  •*•  =  0 
where: 

AVO  =  cascade  m's  violation  of  the  Kutta  condition  for  the 
m 

0  degree  onset  flow. 

AV90  =  cascade  m's  violation  of  the  Kutta  condition  for  the 
m 

90  degree  onset  flow. 

AV  =  cascade  m's  violation  of  the  Kutta  condition  for  the 
m,n 

onset  flow  created  by  the  circulation  on  cascade  n. 

NC  =  number  of  cascades  with  circulation. 
a  =  cascade  average  flow  angle  =  tan 

If  a  is  the  given  cascade  operating  parameter.  Equation  (14)  can  be 

solved  for  the  unknown  values  of  T  .  For  this  case.  Equation  (13)  is 

m 

rewritten  in  matrix  form.  Equation  (14),  which  can  then  be  solved  by 
matrix  reduction.  For  the  usual  case  of  one  cascade,  the  solution  is 


tan  a,  +  tan  cl, 
I  E 


trivial . 
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Avll 

AV12 

avi,nc 

rl 

AVO^  cos  a  +  AV90^  sin  a 

AV21 

AV22  • • * 

AV2,NC 

r2 

AVO2  cos  a  +  AV90 2  sin  a 

• 

• 

• 

• 

=  -u 

avnc,i 

AVNC,2  ”• 

avnc,nc 

rNC 

AVO  cos  a  +  AV90  sin  a 

IN  IN 

For  the  cases  when  a  is  not  the  given  operating  parameter,  Equation  (14) 

cannot  be  used  directly  tc  solve  for  the  T  's;  because,  explicit  relations 

m 


(14) 


do  not  exist  between  a  and  a ,  C.  ,  or  Aa.  When  a_  is  the  given  input 
parameter,  a  matrix  equation  which  contains  the  additional  unknown  a  is 
solved.  For  the  other  cases  when  or  Aa  are  given,  the  program  must 
iterate  to  obtain  the  additional  unknown  a. 

To  obtain  the  matrix  equation  which  contains  a  as  an  unknown  and 
a^.  as  a  known  parameter,  an  additional  equation  is  required.  By  examining 
the  inlet  velocity  triangle  in  Figure  2,  the  following  equation  can  be 
obtained : 


U  cos  a  tan  aT  -  U  sin  a  =  V  (15) 

I  up 

The  upwash  velocity,  V  ,  in  Equation  (15)  is  unknown.  However,  since 

V  is  the  induced  velocity  due  to  the  circulation  on  the  airfoils,  it 
up 

is  known  that  £rm  =  2SP  V  .  The  upwash  velocity  can  therefore  be 
eliminated  from  Equation  (15)  and  after  rearranging  it  yields  the  following: 

tan  a  +  [T^/U  cos  a  +  F^/U  cos  a  ...  r^A1  cos  a]/2SP  =  tan  (16) 

Modifying  Equation  (14)  algebraically  and  inserting  Equation  (16) 
provides  the  necessary  matrix  equation  which  can  be  solved  for  a  and  the 
f  's  for  a  given  a  . 
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UAV901 

Avn 

AV12  • • • 

avi,nc 

UAV902 

AV21 

AV22  * "  • 

AV2,NC 

uav90nc 

%C,1 

AVNC,2 

Avnc,nc 

1 

1 

1 

X 

2SP 

2SP 

2SP 

—  — 

tan  a 

-UAV01 

r  /U  cos  a 

-uavo2 

= 

• 

rNC_i/U  cos  a 

-uavonc 

rNC/U  cos  a 

tan 

—  _ 

_ _  ... 

(17) 


If  Act  or  are  given  as  the  input  parameters,  the  program  will 
iterate  using  an  assumed  value  of  Ct  in  Equation  (14),  until  the  a  is 
found  which  yields  the  given  Aa  or  C  .  Therefore,  equations  which  relate 

Lj 

and  Aa  to  the  I^'s  found  for  a  particular  a  are  required.  These 
equations  are: 


C 


L 


2_ 

Uc 


NC 

l 


m=l 


r 

m 


(18) 


and 


tan  Aa  = 


(19) 


Equation  (18)  is  derived  from  the  definition  of  CL  and  the  basic 

i-j 

relationship  between  circulation  and  the  lifting  force.  Equation  (19) 

can  be  derived  by  writing  tan  Aa  as  tan(a  -  a  ),  expanding  tan(a  -  a  ) 

-L  1  t 

into  an  equation  with  tan  a^  and  tan  a^,  and  then  using  the  velocity 

triangle  to  express  tan  aT  and  tan  a_  in  terms  of  F  ,  sin  a  and  cos  a. 

It  m 
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2 . 3  Transformation  of  Circular  Cascades 

Up  to  this  point,  the  analysis  has  dealt  with  linear  infinite  cascades, 
and  is  applicable  to  the  blade  sections  of  axia]  flow  machines.  For  radial 
flow  machines,  it  is  also  desirable  to  be  able  to  analyze  inlet  guide 
vanes  of  turbines  and  diffuser  vanes  of  pumps  and  compressors.  In  the 
preceding  analysis,  linear  cascades  are  handled  by  the  concept  that 
Equations  (2)  and  (3)  were  derived  for  the  velocities  induced  by  infinite 
arrays  of  sources  and  vortices  of  constant  spacing  SP.  The  flow  field 
thus  repeats  every  spacing  SP  and  satisfying  boundary  conditions  on  one 
airfoil  of  a  cascade  automatically  satisfies  the  boundary  conditions  on 
all  the  airfoils  of  that  cascade. 

A  circular  cascade  can  also  be  considered  an  infinite  cascade; 
because  we  can  continually  observe  points  around  the  cascade,  and  if 
all  the  vanes  are  identical  the  flow  pattern  repeats  after  each  vane. 

If  all  the  vanes  are  different,  the  vanes  can  be  considered  as  several 
infinite  cascades  with  the  flow  pattern  repeating  after  each  revolution. 

Rather  than  deriving  equations  similar  to  Equations  (2)  and  (3)  for 
circular  cascades,  the  approach  taken  was  to  mathematically  transform 
the  circular  cascades  into  linear  cascades  and  then  use  the  existing 
analysis.  Since  the  flow  is  being  analyzed  using  potential  theory,  a 
conformal  transformation  can  be  used  and  the  same  governing  equations 
and  boundary  conditions  will  apply  in  both  the  real  and  transformed 
planes.  The  transformation  used,  is  performed  by  the  following  equations 
and  is  illustrated  in  Figure  3: 

y  =  0  (20) 


x 


In  r 


(21) 
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In  addition  to  the  equations  for  transforming  the  circular  cascade 
geometry  into  a  linear  cascade,  relationships  are  also  required  for 
relating  the  velocities  in  the  two  planes.  From  the  velocities,  the 
pressure  field  can  then  be  calculated.  The  necessary  relationships  can 
be  derived  by  first  defining  a  stream  function  0  in  the  r-9  plane  such 
that: 


and 


v  -iii 

r  r  3e 


V 


e 


30 

3r 


(22) 


(23) 


It  is  also  necessary  to  define  the  stream  function  in  the  x-y  plane  such 
that: 


V 

x 


3  Q 

3y 


(24) 


V  =  -  ^ 


30 

3x 


(25) 


From  the  chain  rule  we  know  that 


and 


3Q  _  30_  _3v 
39  "  3y  38 


3Q  _  d\fr  _3x 
3r  3x  3r 


(26) 


(27) 


Differentiating  Equations  (20)  and  (21)  we  find  that  =  1  and 

o  0 

X  1 

=  — .  Substituting  these  relations  and  Equations  (22)  through  (25) 
into  Equations  (26)  and  (27)  yields  the  desired  relationship  between 
the  velocity  components  in  the  real  and  transformed  planes. 


V  =  -  V 
r  r  x 


(28) 


V0  =  7  Vv 


(29) 
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By  combining  the  velocity  components  into  the  resultant  velocity  a  similar 
equation  relating  the  resultant  velocities  is  obtained: 


re 


•\fl  '  ve  -vf 


2  12 

v  )  +  (-  v  y 

x  r  y 


2  2 

V  +  V 
x  y 


=  -  V 
r  xy 


(30) 


III.  PROGRAM  DESCRIPTION 


The  purpose  of  this  section  is  to  describe  the  key  parameters  and 
variables  found  in  the  program  input  and  output,  discuss  the  types  of 
cascades  which  can  be  analyzed,  and  explain  the  structure  of  the  program 
itself.  The  input  format  required  to  run  the  program  is  not  presented 
in  this  section,  but  is  presented  in  the  users  guide  in  Appendix  A. 

In  part  one  of  this  section,  the  program  parameters  and  variables  are 
defined  and  their  application  to  the  different  types  of  cascades  are 
discussed.  The  parameters  used  to  nondimensionalize  the  data  are  also 
described.  In  part  two,  the  restrictions  which  apply  to  the  three  basic 
types  of  cascades  are  presented,  and  ways  for  circumventing  these 
restrictions  to  enable  some  unusual  cascade  geometries  to  be  analyzed 
are  explained.  In  the  final  part  of  this  section,  the  structure  of  the 
program  and  the  functions  of  each  of  the  major  subroutines  are  described. 
3 . 1  Definition  of  Key  Input  and  Output  Variables 

A  convenient  place  to  start  the  discussion  of  input  and  output 
variables  is  with  the  angles  which  describe  the  flow  and  the  orienta¬ 
tion  of  the  airfoils  in  the  cascades.  The  flow  angles  and  the  stagger 
angle  £  for  a  linear  cascade  are  illustrated  in  Figure  2.  In  terms 
of  the  velocity  components,  the  inlet  angle  and  the  exit  angle 


can  be  defined  as: 
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a 


I 


tan 


-1 


(31) 


aE  =  tan 


-1 


r  ys 


v 


(32) 


For  a  circular  cascade,  the  definitions  are  similar  with  the  x  and  y 
components  of  velocity  replaced  by  the  r  and  6  components,  respectively. 
Thus,  for  a  circular  cascade  with  the  radial  flow  direction  inward  or 


outward,  the  inlet  and  exit  flow  angles  are  defined  as: 

V 

-1  f  9I^ 

VJ 


a  =  tan 


(33) 


a£  =  tan 


-1 


r1] 

r  ' 


(34) 


Equations  (31)  through  (34)  yield  the  sign  which  is  consistent  with  the 
chosen  convention.  All  flow  angles  for  both  linear  and  circular  cascades 
are  positive  in  the  counterclockwise  direction.  The  stagger  angle  for 
circular  cascades  is  also  measured  in  a  counterclockwise  direction,  while 
the  stagger  angle  for  linear  cascades  is  the  only  angle  measured  in  a 
clockwise  direction.  The  opposite  sign  convention  was  chosen  for  the 
linear  cascade  stagger  angle  to  be  consistent  with  the  original  Douglas 
Neumann  program.  For  a  typical  cascade  configuration  as  in  Figure  2, 
the  chosen  convention  also  yields  a  positive  lift  corresponding  to  a 
positive  stagger  angle  when  the  cascade  is  subjected  to  an  axial  inlet 
flow. 

Two  other  angles  which  are  sometimes  used  as  input  to  the  Douglas 
Neumann  program  are  the  flow  turning  angle  Act  and  the  cascade  mean  flow 
angle  a.  These  angles  are  defined  as  follows: 


-28- 


26  June  1981 
AMY :  cag 


Aoi  =  ax  -  a£  (35) 

1  f tan  a  +  tan  aj,\ 

a  =  tan"1  - - - “  (36) 

Tan  a  and  tan  a  in  Equation  (36)  can  be  expressed  in  terms  of  the  x  and 
X  Lj 

y  velocity  components  for  a  linear  cascade.  Making  these  substitutions 

would  show  that  a  is  the  angle  of  the  average  velocity  vector  through  the 

cascades,  and  thus  a  has  physical  significance  for  linear  cascades.  For 

circular  cascades,  the  physical  significance  of  a  is  not  as  clear.  Since 

the  magnitudes  of  the  inlet  and  exit  velocities  for  circular  cascades 

depend  on  the  radius,  no  unique  average  velocity  can  be  defined  without 

specifying  the  radii.  The  inlet  and  exit  velocities  can  each  be  converted 

to  equivalent  velocities  at  a  common  radius  by  applying  conservation  of 

mass  (rV  =  constant)  and  conservation  of  angular  momentum  (rV.  =  constant), 
r  v 

The  angle  of  the  average  of  these  equivalent  velocities  is  a.  At  the 

common  radius  Vr  is  equal  to  Vr  ,  just  as  Vx  is  equal  to  VY  for  a 
I  E  AI  XE 

linear  cascade. 

The  angles  a,  a  ,  a  and  Act  describe  the  overall  characteristics 

X  hi 

of  the  flow  through  the  cascade.  However,  the  bulk  of  the  output  from 
the  program  is  the  velocity  and  pressure  on  the  airfoil  surfaces  and  at 
specified  off-body  points  in  the  flow  field.  The  velocities  and  pressures 
are  both  presented  in  nondimensional  form,  thus,  the  normalizing  velocity 
for  the  various  types  of  cascades  should  be  discussed.  For  linear  cascades, 
the  magnitude  U  of  the  average  velocity  vector  V  is  used  for  non- 
dimensionalizing  velocity  and  pressure.  The  vector  Vg  is  illustrated 
in  Figure  2.  From  the  figure  it  can  be  seen  that: 


- 
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cos  a. 


cos  a 


=  V 


I  cos  a 


=  V 


COS  Olj 
E  cos  a 


(37) 


For  a  circular  cascade,  an  average  velocity  vector  is  difficult  to 
define,  therefore,  a  normalizing  velocity  with  more  physical  significance 
was  selected.  With  uniform  flow,  the  volume  flow  rate  Q  for  a  radial 
flow  system  is  V  2ir  rW.  Even  though  the  flow  is  not  uniform  through  a 
cascade,  the  relationship  for  Q  can  be  used  to  define  a  meaningful  reference 
velocity  by  first  writing: 


V^r  =  Q/2ttW  =  constant  (38) 

Since  V^r  is  equal  to  a  constant,  selecting  a  reference  radius  provides 
a  reference  velocity  which  can  be  used  to  nondimensionalize  the  data  for 
circular  cascades.  The  quantity,  W,  is  the  width  of  the  radial  channel. 


1  Q  _  constant 
r 


r  r  2ttW 
ref 


’ref 


(39) 


The  value  r  ^  is  given  as  input  to  the  program  and  was  selected  as  the 
radius  which  positions  the  airfoils  in  a  cascade.  The  users  guide  in 
Appendix  A  will  further  illustrate  r 

Equation  (39)  defines  the  reference  velocity  for  a  circular  cascade, 
but  it  is  recalled  that  the  flow  field  solution  is  originally  obtained 
in  a  transformed  plane  where  the  reference  velocity  is  U.  A  relationship 

between  U  in  the  transformed  plane  and  Vr  is  therefore  required  in  the 

ref 


program. 


Substituting  r  ^ 


for  r  and  V 

rref 


for  V^.  in  Equation  (28)  yields: 


(40) 


■ 
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From  Equation  (37)  it  is  known  that  Vx  =  U  cos  a,  thus: 

V  =  — —  U  cos  a  (41) 

rref  rref 


Equation  (30),  which  transforms  a  velocity  in  the  x-y  plane  to  a  velocity 
in  the  r-0  plane,  can  now  be  written  in  nondimens ional  form  as: 


r  ,  V 
ref  xy 

r  U  cos  a 


(42) 


Equations  similar  to  Equation  (42)  can  be  written  for  the  velocity 

components  V  and  V  by  dividing  Equations  (28)  and  (29)  by  Equation  (41) . 
r  t) 

All  velocities  in  the  output  of  the  current  version  of  the  Douglas 

Neumann  program  are  nondimensionalized  by  V  and  U  for  circular  or 

rref 

linear  cascades,  respectively.  The  pressures  on  the  body  surfaces  and 

at  off-body  points  are  also  nondimensionalized  by  V  and  U.  The 

rref 


pressure  coefficient  is  defined  as: 


C 

P 


1  - 


ref 


s  ref 


V 


ref 


(43) 


where  V  ,  =  U  for  linear  cascades 
ref 

V  ,  =  V  for  circular  cascades 
ref  r  , 
ref 

As  seen  from  Equation  (43),  the  pressure  coefficient  can  be  calculated 

from  the  local  velocity  V  at  the  point  under  consideration,  but  C  also 

Z  P 

represents  the  difference  between  the  local  static  pressure  and  a 

reference  static  pressure  P  The  reference  static  pressure  is 

s  ref 

defined  as  follows: 


1  2 

Ps  ref  =  PT  "  2  P  Vref 


(44) 
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where  V  ,.  =  U  or  V 

ref  r  .. 

ref 

P  =  total  pressure,  which  is  constant  throughout  the  potential 
flow  field 

The  remaining  parameters  to  be  defined  again  involve  overall  cascade 

properties  rather  than  local  pressures  or  velocities.  For  linear  cascades, 

a  lift  coefficient  C,  x  and  y  force  coefficients  C  and  C  ,  and  a 

L  r  F 

x  y 

moment  coefficient  C,  need  to  be  defined.  The  lift  coefficient  CL 

M  t, 

represents  the  total  lift  of  all  the  cascades;  and  if  it  is  not  an 

input  parameter,  the  program  calculates  from  the  total  circulation 

of  the  one  or  more  cascades.  If  only  one  cascade  is  being  analyzed, 

is  printed  at  the  beginning  of  each  page  along  with  the  other 

coefficients.  However,  for  a  multi-cascade  analysis,  is  printed  at 

the  end  of  the  output  of  all  the  cascades  so  that  it  is  not  interpreted 

as  the  lift  coefficient  of  a  particular  cascade. 

The  circulation  on  a  cascade  induces  an  upwash  and  downwash  on 

the  flow  relative  to  the  average  flow  through  the  cascade.  For  this 

reason,  the  lift  force  is  perpendicular  to  the  average  velocity  vector 

V  .  CT  is  defined  by  the  following  equation  and  is  calculated  in  the 
avg  L 

program  by  the  subsequent  equation: 


C 


L 


■|pU2cW 


(45) 


where  F^  =  sum  of  the  lift  forces  for  one  airfoil  of  each  cascade. 

F.  is  perpendicular  to  V 
L  r  avg 

c  =  chord  =  1 


W  =  span,  (W  =  1  so  that  CL  represents  the  lift  per  unit  span) 

Ld 


A 
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C 


L 


2_ 

Uc 


NC 

1 


m=l 


r 

m 


(46) 


where  T  =  circulation  on  an  airfoil  in  cascade  m 
m 

NC  =  number  of  cascades 

c  =  chord  =  1 

The  coefficients  C^,  ,  C^,  and  are  obtained  in  a  different  manner 
x  y 

than  C^,  because  they  are  computed  by  integrating  the  pressure  distributions 

around  the  airfoils.  For  a  multicascade  analysis,  coefficients  are  computed 

for  each  cascade.  C  and  C„  nondimens ionally  represent  the  forces  on 

r  F 

x  y 

the  airfoils  in  the  x  and  y  directions,  respectively,  in  the  cascade  global 

coordinate  system.  The  moment  represented  by  C  is  computed  about  a  point 

M 

specified  in  the  input.  These  coefficients  are  defined  as: 


F 

x 


1  2 

f>trCW 


(47) 


F 

_ 2L_ 

|pU2cW 


(48) 


M 

^  "  ^TT2  2TT 

-^pU  c  W 


(49) 


where  Fx  =  force  on  the  airfoil  in  the  global  x  direction 

Fy  =  force  on  the  airfoil  in  the  global  y  direction 

M  =  moment  about  a  point  specified  in  the  input  (a  positive 

moment  is  counterclockwise) 

c  =  chord  of  the  airfoils  in  the  cascade  for  which  the 


coefficients  are  calculated 


pi 


m 
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C  is  related  to  C_  and  C  through  the  mean  flow  angle  a.  If 
L  r  r 

x  y 

there  is  more  than  one  cascade,  the  force  coefficients  can  be  summed 

and  used  to  compute  C  . 

Li 

NC  NC 

C  =  £  Cp  c  cos  a  -  £  Cp  c  sin  a  (50) 

m=l  m=l  xm 


Equation  (50)  is  not  used  in  the  program,  but  it  is  provided  to  explain 

the  relationship  between  Cp^,  Cp^  and  C^.  If  is  hand  calculated  from 

Equation  (50),  the  resulting  value  should  only  differ  from  the  computer 

value  by  the  small  error  incurred  when  integrating  the  pressure  distributions 

for  Cp  and  Cp  .  The  chord  c  is  required  in  Equation  (50)  for  multiple 
x  y 

cascade  analyses  where  all  the  airfoils  may  not  have  c=l. 

For  circular  cascades,  is  not  computed  because  an  average 
velocity  vector  is  not  easily  defined  and  would  have  little  meaning. 

Two  force  coefficients  and  a  moment  coefficient  are  computed  for  each 
circular  cascade  by  integrating  the  pressure  distributions.  The  following 
equations  define  these  coefficients  and  complete  the  definition  and 
discussion  of  the  input  and  output  variables. 


C 


M 


F 

r 


1  „2 

7  pVv 

cW 

ref 

Fe 

cW 

2  rref 

m 

~  pV2  c^ 

2  rref 


(51) 


(52) 


(53) 
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where  F  =  force  on  the  airfoil  in  the  r  direction 
r 

F  =  force  on  the  airfoil  in  the  9  direction 

y 

M  =  moment  about  a  point  specified  in  the  input  (a  positive 
moment  is  counterclockwise 
3. 2  Cascade  Modeling 

There  are  three  basic  types  of  cascades  which  the  current  version  of 
the  Douglas  Neumann  cascade  program  can  handle.  The  program  can  analyze 
linear  cascades,  circular  cascades  with  the  radial  flow  inward,  and 
circular  cascades  with  the  radial  flow  outward.  In  addition  to  single 
cascade  analyses,  multiple  cascades  can  be  analyzed  as  a  system.  For 
each  of  the  various  types  of  cascade  analyses,  the  following  restrictions 
apply  to  the  cascades: 

1.  All  the  bodies  of  a  particular  cascade  must  be  equally  spaced. 

2.  All  the  bodies  of  a  particular  cascade  must  be  identical. 

3.  For  an  analysis  with  multiple  cascades,  all  the  cascades  must 
have  the  same  spacing  between  the  bodies. 

4.  No  cascade  can  move  relative  to  another  cascade. 

The  last  restriction  is  the  only  restriction  which  cannot  be  overcome 
through  cascade  modeling.  Under  no  circumstance  can  a  cascade  move 
relative  to  another  cascade  because  this  would  constitute  an  unsteady 
flow  field. 

For  simple  cascade  geometries,  no  modeling  is  involved  in  the  program 
input.  The  user  simply  inputs  to  the  program  the  airfoil  geometry  and 
cascade  parameters  such  as  spacing,  stagger  angle,  inlet  flow  angle, 
etc.  However,  for  unusual  cascade  geometries  where  one  of  the  first 
three  restrictions  is  violated,  a  cascade  can  often  be  modeled  as  several 
cascades  to  overcome  the  restrictions. 


] 
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Two  examples  will  be  used  to  illustrate  cascade  modeling.  For  the 
first  example,  the  cascade  of  circular  cylinders  shown  in  Figure  4  has 
unequal  spacing  between  the  cylinders  and  thus  violates  the  first 
restriction.  In  this  case,  although  the  distances  from  one  cylinder 
to  the  next  cylinders  above  and  below  are  not  equal,  it  is  seen  that  the 
cascade  pattern  repeats  every  spacing  SP.  The  single  cascade  of  unequally 
spaced  cylinders  can  thus  be  modeled  as  two  cascades.  Both  cascades  will 
have  the  same  spacing  between  the  cylinders,  with  the  second  cascade 
offset  from  the  first  cascade  by  the  distance  ADDY  shown  in  the  figure. 

The  numbers  on  the  cylinders  in  Figure  4  identify  the  cylinders  belonging 
to  the  two  cascades. 

Although  in  the  first  example  the  cylinders  of  cascade  one  and  two 
have  the  same  diameter,  the  two  cascades  could  have  different  size  or 
differently  shaped  bodies.  Thus,  the  same  procedure  can  be  used  to 
model  cascades  with  bodies  of  different  geometries.  The  process  of 
representing  differently  shaped  bodies  or  equally  spaced  bodies  as  different 
cascades  could  be  carried  to  an  extreme  where  a  random  cluster  of  differently 
shaped  two-dimensional  objects  could  be  analyzed.  Each  object  could  be 
represented  as  a  cascade,  and  the  spacing  could  be  made  very  large  such 
that  the  effect  of  other  members  of  the  same  cascade  would  be  negligible. 

This  extreme  situation  may  require  a  very  large  number  of  points  to  represent 
the  different  body  geometries  and,  thus,  may  not  be  practical.  However, 
the  extreme  case  does  illustrate  the  flexibility  of  the  program  made 
available  by  the  ability  of  the  program  to  analyze  multiple  cascades. 

It  should  be  possible  to  model  most  cascades  of  engineering  significance 
with  a  reasonably  number  of  cascades  which  will  not  violate  any  of  the 


program  restrictions. 
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The  second  example  of  cascade  modeling  illustrates  an  analysis  of 
a  cascade  with  apparently  differently  shaped  airfoils.  Figure  5 
is  a  schematic  of  stay  vanes  and  wicket  gates  of  a  hydraulic  turbine. 

It  is  noticed  from  the  figure  that  the  stay  vanes  have  five  different 
geometries.  Over  most  of  the  circumference  the  stay  vane  geometry 
is  constant  with  either  stay  vane  A  or  stay  vane  B  being  employed.  In 
the  regions  where  the  stay  vane  geometry  does  not  vary,  the  flow  can 
easily  be  analyzed  with  the  current  version  of  the  Douglas  Neumann 
program.  The  wicket  geometry  and  the  appropriate  stay  vane  geometry 
are  input  to  the  program,  neglecting  the  other  stay  vane  geometries 
which  are  remote  from  the  area  of  interest. 

Neglecting  the  variation  of  the  stay  vane  geometry,  the  wicket 
gates  and  stay  vanes  represent  two  cascades  which  do  not  violate  any 
of  the  program  restrictions.  Suppose,  however,  that  it  is  desirable 
to  analyze  the  flow  in  the  region  of  stay  vane  D,  where  it  is  seen 
from  Figure  5  that  there  are  five  different  stay  vane  geometries  side 
by  side.  Since  there  is  no  symmetry  in  the  distribution  of  stay  vane 
geometries,  to  analyze  the  entire  flow  field  including  the  five  different 
stay  vane  geometries  would  require  48  cascades.  Each  wicket  gate  and 
stay  vane  would  need  to  be  represented  as  a  separate  cascade  with  the 
spacing  SP=360° .  The  number  of  points  required  to  represent  48  cascades 
certainly  makes  this  approach  unrealistic. 

A  more  realistic  approach  to  analyze  the  flow  near  vane  D  is  to 
only  consider  the  two  vanes  on  both  sides  of  vane  D,  and  after  these 
vanes  allow  the  pattern  to  repeat.  In  this  way,  the  effects  of  the  two 
nearest  vanes  in  both  directions  is  included,  but  the  effects  of  the  true 


vane  shapes  beyond  the  first  two  vanes  are  neglected  in  the  model.  This 
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approach  is  illustrated  in  Figure  6.  In  order  that  the  wicket  gate 
spacing  be  the  same  as  the  stay  vane  spacing,  the  wicket  gates  corresponding 
to  the  five  stay  vanes  must  also  be  modeled  as  different  cascades.  The 
need  to  represent  the  wicket  gates  as  five  cascades,  even  though  they 
have  the  same  geometry,  is  also  apparent  from  the  fact  that  they  are 
each  downstream  of  a  different  stay  vane  geometry  and  will  each  have  a 
different  flow  solution.  The  flow  in  the  region  of  stay  vane  D  can 
therefore  be  modeled  with  ten  cascades,  the  maximum  number  the  program  is 
currently  set  up  to  handle. 

In  Figure  6  the  ten  cascades  representing  the  wicket  gates  and  stay 
vanes  near  vane  D  are  numbered  to  identify  the  members  of  the  various 
cascades.  The  original  vane  shapes  are  also  identified  by  the  letters. 

It  is  seen  that  each  cascade  has  the  same  spacing  SP  and  after  the  spacing 
SP  the  pattern  repeats.  The  relative  position  of  the  various  cascades  is 
specified  by  the  radius  and  the  angles  ADDY^,  ADDY^,  etc. 

It  is  important  to  remember  when  modeling  different  vane  shapes, 
as  in  the  preceding  example,  that  the  vane  of  primary  concern  should  be 
kept  in  the  middle  of  the  group.  This  will  minimize  errors  caused  by 
not  having  the  true  vane  shape  a  few  vanes  from  the  place  of  concern. 

For  example,  the  model  shown  in  Figure  6  should  not  be  used  to  analyze 
the  flow  near  vane  C,  but  instead,  a  new  model  with  two  true  vanes  on 
both  sides  of  vane  C  should  be  employed. 

Note:  A  circular  cascade  becomes  an  infinite  linear  cascade  in  the 

transformed  plane  because  each  succeeding  revolution  around  the  circular 
cascade  will  generate  more  of  the  same  airfoils.  In  the  example  where 
there  were  24  stay  vanes,  the  true  pattern  would  repeat  every  24  vanes. 


However,  the  stay  vanes  were  modeled  as  5  cascades  which  cannot  be  evenly 
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Figure  6 


Illustration  of  One  Possible  Wav  to  Model  the  Wicket 
Oates  and  Stay  Vanes  in  the  Region  of  Stay  Vane  D 
of  Figure  5 

i  ^ MM— I 


divided  into  24.  Therefore,  the  pattern  will  not  repeat  every  24  vanes 
as  it  should.  This  fact  may  be  somewhat  disturbing  because  transforming 
the  linear  cascade  back  to  the  circular  cascade  will  result  in  several 
different  vane  shapes  occupying  the  same  location.  In  the  transformed 
plane,  where  the  flow  field  solution  is  originally  obtained,  the  airfoils 
which  will  occupy  the  same  location  in  the  real  plane  are  24  vanes  apart. 
Thus,  they  have  essentially  no  effect  on  each  other,  and  the  fact  that 
the  pattern  does  not  repeat  every  24  vanes  is  inconsequential.  The 
author  has  not  verified  this  conclusion  with  computer  runs  of  the  program, 
therefore,  it  may  be  wise  for  the  user  to  do  so  if  the  number  of  cascades 
used  in  modeling  a  circular  cascade  cannot  be  divided  evenly  into  the 
total  number  of  vanes. 

3.3  Structure  of  the  Program 

The  current  version  of  the  Douglas  Neumann  cascade  program  consists 
of  a  main  program  and  five  major  subroutines.  The  main  program's  only 
function  is  to  call  the  appropriate  major  subroutines.  Each  of  the  major 
subroutines  will  be  described  in  the  following  paragraphs.  Besides  the 
CALL  statements  for  the  subroutines,  the  main  program  also  contains  two 
logical  IF  statements.  One  IF  statement  selects  SUBROUTINE  PART3  or 
SUBROUTINE  PART4  for  solving  the  set  of  simultaneous  equations  depending 
on  the  size  of  the  problem.  The  other  IF  statement  directs  the  program 
to  go  to  SUBROUTINE  PART5  following  SUBROUTINE  PARTI  for  cases  where  the 
basic  solutions  of  the  previous  case  can  be  used  for  the  present  case.  The 
main  program  will  continue  to  process  cases  until  an  input  code  directs  the 
program  to  be  terminated.  If  cases  of  similar  geometry  are  run  where 
only  the  flow  angle  or  lift  coefficient  is  changed,  the  same  basic 
solutions  can  be  used  for  each  case.  Subsequent  to  the  first  solution 
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when  the  basic  solutions  are  obtained,  the  program  only  needs  to  solve 
the  combination  equations.  This  feature  will  save  a  considerable  amount 
of  computation  time  when  several  solutions  with  the  same  cascade  geometry 
but  different  flow  conditions  are  required: 

SUBROUTINE  PARTI: 

The  major  purpose  of  SUBROUTINE  PARTI  is  to  read  the  input  data 
and  do  any  of  the  preliminary  calculations  which  are  necessary  before 
the  set  of  simultaneous  equations  can  be  set  up  and  solved.  One  of 
the  first  cards  to  be  read  is  a  control  card  specifying  the  type  of 
cascade  and  giving  other  run  control  information.  Cascade  parameters 
are  read  next,  followed  by  the  airfoil  coordinate  data.  The  airfoil 
coordinate  data  is  usually  read  from  cards.  However,  if  the  airfoil 
geometry  is  the  same  as  the  previous  case,  the  data  can  be  read  from 
a  file.  For  runs  where  a  particular  airfoil  geometry  is  repeated  for 
modeling  purposes,  as  in  the  stay  vane  -  wicket  gate  example  in  the  preceding 
subsection,  SUBROUTINE  PARTI  will  generate  the  airfoil  coordinate  data 
for  the  similar  airfoils  following  the  input  of  the  first  airfoil  data. 

After  the  coordinate  data  have  been  input,  the  data  are  stored  in  a  file 
before  the  required  maniuplations  begin. 

For  circular  cascades,  SUBROUTINE  PARTI  reorders  the  coordinate 
data  so  that  in  the  transformed  plane  the  data  will  have  the  same  order 
as  data  for  linear  cascades.  The  necessity  of  reordering  the  data  is 
described  in  more  detail  in  Appendix  C.  For  circular  cascades  with  the 
radial  flow  direction  inward,  the  coordinate  data  are  put  in  reverse 
order.  If  the  radial  flow  direction  is  outward,  the  starting  point  of 
the  data  is  shifted  to  the  trailing  edge  of  the  airfoil  for  the  current 


flow  direction. 
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SUBROUTINE  PARTI  next  scales  the  coordinate  data  by  the  value  of 
the  chord  given  for  each  cascade.  The  coordinate  data  are  then  converted 
from  local  to  global  coordinates,  which  involves  translating  and  rotating 
the  coordinates  so  the  airfoils  have  the  proper  location  and  orientation 
in  the  flow  field.  For  circular  cascades,  the  transformation  is  carried 
out  at  the  same  time  the  coordinates  are  converted  to  global  coordinates. 
Coordinates  of  the  element  midpoints  are  determined  and  the  sine  and  cosine 
of  the  angle  of  the  elements  are  calculated.  The  airfoil  coordinates  are 
put  in  complex  form  to  be  used  in  the  complex  arithmetic  involved  in 
setting  up  the  A  matrix  in  Equation  (11). 

The  last  section  of  SUBROUTINE  PARTI  reads  in  or  generates  the 
coordinates  of  off-body  points  where  the  flow  is  to  be  determined.  The 
coordinates  of  each  off-body  point  can  simply  be  read  in  from  a  card. 
However,  in  many  cases  it  is  desirable  to  locate  a  prescribed  number  of 
off-body  points  equally  spaced  across  a  section  of  the  flow  at  a  constant 
x  or  r  value.  Under  these  circumstances,  the  program  can  generate  the 
coordinates  of  the  off-body  points  after  being  given  the  initial  data 
for  a  set  of  points.  By  supplying  the  program  with  the  coordinates  of 
the  first  point  of  the  set,  the  number  of  off-body  points  in  the  set, 
and  the  increment  Ay  or  A8  by  which  to  locate  the  remaining  points,  the 
program  will  generate  the  required  coordinates.  Any  number  of  sets  of 
off-body  points  can  be  used,  provided  the  total  number  of  off-body  points 

does  not  exceed  the  array  dimensions.  For  circular  cascades,  the  off-body 
points  are  transformed  to  the  linear  cascade  by  SUBROUTINE  PARTI. 
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SUBROUTINE  PART 2: 

In  SUBROUTINE  PART2,  the  set  of  simultaneous  equations  which  are 

solved  for  the  source  distributions  on  the  airfoils  are  set  up.  Setting 

up  the  equations  consists  of  calculating  the  values  of  the  A  matrix 

and  determining  the  onset  flows  in  Equation  (11) .  Complex  arithmetic 

is  used  to  calculate  the  A  matrix,  such  that  the  real  and  imaginary  parts 

yield  the  coefficients  for  the  flow  induced  both  tangent  and  normal  to 

the  airfoil  surface.  The  coefficients  representing  the  flow  induced  normal 

to  the  surface  are  the  values  in  the  A  matrix.  The  coefficients  for 

the  flow  induced  tangent  to  the  surface  are  stored  and  later  used  to 

calculate  the  velocity  on  the  airfoil  surfaces  after  the  source  strengths 

are  known.  The  sum  of  the  real  parts  and  the  sum  of  the  imaginary  parts 

of  the  coefficients  for  each  row  of  the  matrix  are  also  calculated  in 

n? 

SUBROUTINE  PART2.  The  sum  of  the  real  parts  is  the  term  £  A'.,  in 

k=nx 

Equation  (12)  and  provides  the  onset  flows  induced  normal  to  the  element 
surfaces  by  circulation.  The  sums  of  the  imaginary  parts  provide  the 
velocities  induced  tangent  to  the  elements  by  circulation  and  are  required 
to  calculate  the  final  velocities  at  the  element  midpoints. 

SUBROUTINE  PART 3  and  SUBROUTINE  PART4 : 


The  set  of  simultaneous  equations  for  each  basic  solution  are  solved 
in  either  SUBROUTINE  PART3  or  SUBROUTINE  PART4  depending  on  the  size  of  the 
problem.  SUBROUTINE  PART4  reads  the  coefficients  and  onsets  flows  cal¬ 
culated  in  SUBROUTINE  PART2  from  files  and  stores  the  values  in  arravs. 
SUBROUTINE  MIS1  is  then  called,  which  solves  the  equations  by  matrix 
reduction.  Since  SUBROUTINE  PART4  and  MIS1  require  the  entire  A  matrix 
and  onset  flow  matrices  to  be  in  the  computer  core  at  once,  the  use  of 
these  subroutines  is  restricted  by  the  size  of  the  matrices  which  can  be 


accommodated  in  the  core. 
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There  are  currently  two  versions  of  SUBROUTINE  PART3  which  will 
solve  the  equations  when  it  is  not  practical  or  possible  to  have  the 
entire  A  matrix  in  the  computer  core.  The  original  Douglas  Neumann 
cascade  program  employed  an  iterative  technique  to  solve  the  equations 
for  large  problems.  This  iterative  technique  is  found  in  one  of  the 
two  current  versions  of  SUBROUTINE  PART3.  A  disadvantage  of  the  iterative 
approach  is  that  the  complete  set  of  equations  must  be  read  from  auxiliary 
storage  for  each  iteration.  Reading  data  from  auxiliary  storage  is  one 
of  the  slower  operations  of  computer  processing  and,  thus,  may  greatly 
increase  the  computation  time. 

To  reduce  the  computation  time  for  large  problems,  a  new  solution 
procedure  was  developed  at  ARL/PSU  during  the  current  rewriting  of  the 
Douglas  Neumann  program.  The  new  procedure  found  in  the  second  version 
of  SUBROUTINE  PART3  reduces  the  A  matrix  by  operating  on  blocks  of  rows 
of  the  A  matrix  in  the  core.  The  procedure  is  designed  to  take  advantage 
of  whatever  core  space  is  available  and  it  continually  optimizes  the  equation 
storage  by  increasing  the  number  of  equations  in  the  core  as  the  number  of 
nonzero  terms  in  each  equation  decreases.  Comparative  computer  runs 
solving  the  same  problem  using  the  three  solution  procedures  show  that 
the  new  procedure  requires  the  least  amount  of  computer  time. 

The  three  available  solution  procedures  found  in  SUBROUTINE  PART4 
and  the  two  versions  of  SUBROUTINE  PART3  are  described  in  more  detail 
in  Appendix  B.  Discussions  of  the  storage  requirements  and  the  major 
advantages  and  disadvantages  of  each  procedure  are  included.  Appendix  B 
also  describes  the  program  dimension  statements  so  that  the  user  can 
choose  a  solution  procedure  and  modify  the  array  dimensions  to  match  the 
program  for  a  particular  problem  to  the  computer  system. 


A 
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SUBROUTINE  PART 5: 

The  last  major  subroutine  is  SUBROUTINE  PART5.  Its  first  function 
is  to  set  up  and  solve  the  combination  equations.  The  solution  of  the 
combination  equations  enables  the  basic  solutions  to  be  combined  such  that 
the  cascade  operating  condition  and  the  Kutta  condition  are  satisfied.  As 
discussed  in  Section  2.2,  the  combination  equations  are  based  on  a  AV 
across  the  trailing  edge  which  represents  each  of  the  basic  solution's 
violation  of  the  Kutta  condition.  The  beginning  of  SUBROUTINE  PART5 
calculates  the  required  AV's. 

Following  the  determination  of  the  AV's  there  are  four  sections  of 
the  subroutine  which  solve  the  combination  equations  depending  on  which 
of  the  four  possible  input  parameters  were  given.  If  a  or  a^.  is  the  given 
cascade  operating  parameter,  the  combination  equations  are  solved  by  matrix 
reduction.  If  CL  or  Act  is  given,  the  program  iterates  to  determine  the 
value  of  a  the  flow  must  have  to  yield  the  given  C  or  Act.  The  solution 

Lj 

of  the  combination  equations  produces  the  value  of  the  circulation  for 
each  cascade. 

Following  the  solution  of  the  combination  equations,  SUBROUTINE  PART5 
calculates  the  velocity  at  the  midpoint  of  each  element.  This  is  done  by 
adding  the  onset  velocities  and  the  velocities  induced  by  all  the  elements 
for  each  of  the  basic  solutions.  When  adding  solutions,  the  velocities 
of  the  appropriate  basic  solutions  are  scaled  by  sin  a,  cos  a  or  the  value 
of  the  circulation  for  each  cascade.  For  circular  cascades,  the  velocities 
are  next  transformed  back  to  the  real  plane  and  the  data  are  reordered  to 
correspond  to  the  data  original  input  order.  The  airfoil  pressure  distribu¬ 
tion  is  calculated  from  the  velocities  and  then  is  integrated  to  yield 


-47- 


26  June  1981 
AMY:  cag 


force  and  moment  coefficients.  The  results  are  printed  and  the  velocity 
calculation  process  repeats  for  the  other  cascades  of  a  multi-cascade 
analysis.  The  last  section  of  SUBROUTINE  PARTS  calculates  the  velocity 
components  and  the  pressure  at  each  of  the  specified  off-body  calculation 
points . 

IV.  PROGRAM  VERIFICATION 
4 . 1  Methods  of  Verification 

Every  computer  program  should  be  checked  to  assure  that  it  is  calculating 
what  it  was  intended  to  calculate.  Since  the  Douglas  Neumann  program  yields 
an  approximate  solution  to  a  potential  flow  problem,  the  obvious  means 
of  verifying  the  program  is  to  analyze  the  flow  through  a  cascade  for  which 
an  exact  potential  flow  solution  exists.  In  this  section,  results  from  the 
current  rewritten  version  of  the  Douglas  Neumann  program  are  compared  to 
exact  solutions  for  various  types  of  cascades.  The  exact  solutions  are 
obtained  by  using  conformal  transformations  to  convert  the  potential  flow 
around  a  single  cylinder  to  the  flow  around  airfoils  in  infinite  cascades. 
Besides  the  comparisons  with  exact  solutions,  one  additional  comparison 
is  made  between  the  results  from  the  current  version  of  the  program  and 
results  from  the  original  version  which  were  presented  in  Reference  [1] . 

This  comparison  was  made  to  assure  the  current  version  of  the  program 
can  analyze  multiple  cascades  of  lifting  and  nonlifting  bodies  as  did  the 
original  program. 

It  is  also  important  to  know  how  well  a  computer  program  predicts 


real  phenomena  or  experimental  data.  However,  the  degree  to  which 
calculated  results  match  experimental  data  is  not  necessarily  a  check 
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on  the  proper  functioning  of  a  computer  program.  Instead,  comparisons 
with  experimental  data  are  checks  on  how  well  the  simplifying  assumptions 
and  governing  equations  represent  the  real  situation.  Comparisons  of 
results  from  the  Douglas  Neumann  program  and  experimental  data  have  been 
made  in  other  reports  but  will  not  be  made  here.  In  Reference  [1],  measured 
pressure  distributions  for  linear  cascades  are  compared  to  calculated  results 
from  the  program.  Measured  and  calculated  flow  angles  for  a  circular  cascade 
with  the  radial  flow  direction  inward  are  compared  in  Reference  [2]. 

With  the  current  version  of  the  Douglas  Neumann  program  having  the 
additional  capability  of  analyzing  circular  cascades,  several  new  solutions 
as  well  as  the  original  solutions  for  linear  cascades  required  verifica¬ 
tion.  The  types  of  cascades  for  which  the  solutions  required  verification 
were: 

1.  A  linear  cascade. 

2.  A  circular  cascade  with  the  radial  flow  inward. 

3.  A  circular  cascade  with  the  radial  flow  outward. 

4.  Multiple  cascades  of  each  of  the  first  three  types. 

5.  Multiple  cascades  consisting  of  both  lifting  and  nonlifting  bodies. 
Solutions  for  the  first  three  types  of  cascades  were  verified  by  comparing 
them  with  their  corresponding  exact  analytical  solution.  For  multiple 
cascades  of  each  type,  a  single  cascade  was  modeled  as  several  cascades, 
thus  the  same  exact  solution  could  be  used  to  verify  both  single  and 
multiple  cascades.  As  previously  mentioned,  the  fifth  type  of  cascade 
solution  was  verified  by  comparing  the  results  from  the  current  version 

of  the  program  with  the  results  presented  in  Figure  5  of  Reference  [1]. 
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The  cascade  parameters  (C^  is  a  parameter  for  linear  cascades 
only),  cij,  a,  and  Aa  were  all  known  from  the  exact  solutions  for  the 
various  types  of  cascades.  Thus,  each  of  the  cascade  parameters  could 
be  used  as  input  to  the  program  to  check  the  various  input  options.  Using 
each  of  the  different  cascade  parameters  as  input,  the  other  parameters 
were  calculated  by  the  program.  In  this  way,  the  cascade  parameters 
calculated  by  the  program  were  verified.  As  required,  the  calculated 
pressure  distributions  on  the  airfoils  were  nearly  identical  regardless 
of  the  input  parameter.  Velocity  components  at  points  off  the  bodies 
were  also  obtained  from  the  exact  analytical  solutions  and  were  used  to 
check  the  off-body  data  calculated  by  the  program. 

4 . 2  Exact  Analytical  Cascade  Solutions 

Exact  potential  flow  solutions  for  linear  and  circular  cascades  can 
be  obtained  from  the  flow  around  a  circular  cylinder  by  employing  a  series 
of  conformal  transformations.  The  transformations  used  to  obtain  a  linear 
cascade  are  the  same  as  those  found  in  References  [1]  and  [3].  An  additional 
transformation  was  then  employed  to  transform  the  linear  cascade  to  a 
circular  cascade. 

Flow  around  a  circular  cylinder  is  generated  by  locating  a  pair 
of  sources  and  vortices  at  points  B1  and  S'  and  a  pair  of  sinks  and 
vortices  at  points  A’  and  A"  shown  in  the  S-plane  of  Figure  7.  The 
locations  of  the  source,  sink  and  vortex  pairs  determines  the  position 
of  the  cylinder.  The  position  of  the  cylinder  in  turn  determines  the 
final  shape  of  the  airfoils  after  the  transformation.  To  form  a  cylinder, 
the  sink  and  vortex  at  point  Af  must  mirror  the  sink  and  vortex  at  A", 
and  similarly  the  singularities  at  B'  and  B"  must  do  the  same.  The  two 
points  of  each  pair  lie  on  the  same  radial  line  and  the  radial  distance 


Figure  7.  Illustration  of  the  Transformation  of  a  Cylinder  in  the  S-Plane  - 
Linear  Cascade  in  the  Z-Plane  and  the  Addition  Transformation  of 
the  Circular  Cascade  in  the  T-Plane.  The  Radial  Flow  is  Outward 
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from  the  center  of  the  cylinder  to  the  inner  point  is  the  inverse  of 
the  radial  distance  to  the  outer  point  (i.e.  if  A'  is  r  from  the  center, 
A"  is  1/r  from  the  center).  To  satisfy  continuity,  the  source  and  the 
sink  have  the  same  strength.  Together,  the  source,  sli/k  and  vortex 
strengths  determine  the  inlet  flow  angle  and  must  also  be  specified  so 
as  to  satisfy  the  Kutta  condition.  With  this  brief  background,  the 
following  parameters  and  coordinates  provide  the  necessary  information 
for  the  cylinder  used: 


a  =  radius  of  the  cylinder  =  1.0 

d  =  distance  to  the  point  where  the  cylinder  intersects 

the  xg  axis  (this  point  becomes  the  T.E.  of  the  airfoils) 
d  =  .8460254 


Point 

Center  of  the  Circle 
A' 

A" 

B’ 

B" 


Coord inates 
(-0.02,  0.5) 

(1.06605,  0.0) 
(0.739739,  0.150228) 
(-1.088605,  0.0) 
(-0.787722,  0.140783) 


Two  sequential  transformations  are  used  to  map  the  cylinder  in  the 
complex  S-plane  into  a  cascade  of  airfoils.  The  first  transformation 
maps  the  cylinder  into  a  single  airfoil  shape  in  the  complex  Q-plane 
as  illustrated  in  Figure  7.  This  transfoi-mation  is  performed  by  the 
following  equation: 
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Q  -  rd 
Q  +  rd 


S  -  d 


S  +  d 


(54) 


where  r  =  trailing-edge-angle  constant  =  1.85 
Q  =  complex  coordinates  in  the  Q-plane 
S  =  complex  coordinates  in  the  S-plane 
The  second  transformation  maps  the  airfoil  shape  in  the  complex 
Q-plane  into  a  cascade  of  airfoils  in  the  complex  Z-plane.  This 
transformation  is: 


Z  =  In 


(55) 


A  and  B  are  the  complex  coordinates  of  the  points  in  the  O-plane  c<  rresponding 
to  the  points  A'  and  B'.  Points  A  and  B  are  illustrated  in  Figure  7.  The 
numerical  value  of  A  and  B  can  be  found  by  plugging  the  given  coordinates 
for  A'  and  B'  into  Equation  (54),  respectively. 

The  linear  infinite  cascade  of  airfoils  resulting  from  the  second 
transformation  is  illustrated  in  the  Z-plane  of  Figure  7.  This  is  the 
cascade  which  was  analyzed  with  the  current  version  of  the  Douglas  Neumann 
program  to  verify  the  program  functions  properly.  The  cascade  was  analyzed 
as  a  single  cascade  and  was  also  modeled  and  analyzed  as  multiple  cascades. 

As  an  example  of  a  multiple  cascade  model,  the  numbers  above  the  airfoils 
in  Figure  7  identify  the  members  of  each  cascade  for  a  two  cascade  model. 

To  obtain  a  circular  cascade  for  which  an  exact  potential  flow 
solution  is  possible,  a  third  transformation  is  applied  to  the  linear 
cascade.  Before  the  third  transformation  is  applied,  however,  the  linear 
cascade  is  modified  slightly.  For  the  case  illustrated  in  Figure  7, 
where  the  circular  cascade  has  the  radial  flow  direction  inward,  the 


linear  cascade  is  translated  and  the  signs  of  the  x^,  coordinates 
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are  changed.  These  changes,  resulting  in  the  cascade  illustrated  in 
the  z' -plane  of  Figure  7,  are  necessary  to  obtain  the  desired  airfoil 
orientation  and  flow  direction  for  the  circular  cascade  in  the  T-plane. 
In  Figure  8,  the  transformations  leading  to  a  circular  cascade  with 
the  radial  flow  direction  outward  are  illustrated.  Up  through  the 
Z-plane,  the  transformations  are  identical  to  those  in  Figure  7.  In 
the  Z' -plane  of  Figure  8,  however,  it  is  seen  that  the  cascade  only 
needs  to  be  translated  to  the  +x  direction  to  yield  the  desired  airfoil 
orientation  and  flow  direction  in  the  T-plane. 

The  transformation  used  to  go  from  the  Z' -plane  to  the  T-plane  is 
the  same  transformation  used  in  the  program,  except  that  for  the  exact 
solutions  it  is  used  in  complex  form.  To  obtain  the  new  form.  Equations 
(20)  and  (21)  are  first  repeated  here  in  terms  of  the  nomenclature  of 
Figures  7  and  8. 


V  =  eT 

(56: 

X-  » 

=  In  rT  or  rT  =  e  z 

(57: 

It  is  necessary  to  express  the  cartesian  coordinates  of  the  T-plane 
in  terms  of  r^  and  9^. 


xT  =  rT  cos  0T 


(58: 


yT  =  rT  sin  0T  (59! 

Applying  the  identities  of  Equations  (56)  and  (57)  to  Equations  (58)  and 
(59)  yields: 


Q  PLANE  .  -v  Z  PLANE  ,  *-v  LOCAL  COORDINATE 


Figure  8.  Illustration  of  the  Transformations  Which  Yield  a  Circular  Cascade 
in  the  T-P lane  Witli  the  Radial  Flow  Inward 
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xT  =  eXz'  cos  y  ,  (60) 

yT  =  eXz'  sin  y^,  (61) 

Expressing  x  and  as  a  complex  number  yields  the  new  form  of  the 
transformation. 

T=xT+iyT=ez  (cos  y^,  +  i  sin  y  , ) 

=  eXz'  eiyz '  =  eXz’  +  iyz’ 

Z' 

T  =  e  (62) 

This  new  expression  of  the  transformation  in  complex  form  facilitates 
the  calculation  of  the  complex  velocities  in  the  exact  solution  of  the 
flow  field  for  the  circular  cascades. 

Starting  with  coordinates  of  points  on  the  surface  of  the  cylinder 
in  the  S-plane  the  preceding  transformations  define  the  cascade  geometries 
for  which  an  exact  potential  flow  solution  is  obtained.  The  transformations 
also  provide  the  local  coordinates  of  points  defining  the  airfoil  shapes, 
which  are  required  as  input  to  the  Douglas  Neumann  program.  To  obtain  the 
exact  flow  field  solution,  however,  it  is  necessary  to  apply  the  trans¬ 
formations  to  the  complex  potential  function  which  represents  the  flow 
around  the  cylinder  in  the  S-plane.  The  complex  potential  function  F(S) 
is  derived  by  adding  the  contribution  of  each  singularity 


56’ 


26  June  1981 
AMY : cag 


F(S)  =  -E  ln( S  -  A')  +  iG  ln(S  -  A')  -  E  ln(S  -  A")  -  iG  ln(S  -  A") 

sink  vortex  ^  sink  vortex^ 


+  E  ln( S  -  B’)  -  iH  ln(S  -  B')  +  E  ln(S  -  B")  +  iH  ln(S  -  B") 

source  vortex  source  vortex  ^ 


where  E  =  source  and  sink  strength 

G  =  vortex  strength  at  A'  and  A" 

H  =  vortex  strength  at  B  1  and  B" 

Equation  (63)  can  be  written  in  simplified  form  as: 


F(S)  =  Kx  In  (S  -  A*)  +  K2  In  (S  -  A")  +  K2  In  (S  -  B')  +  K2  In  (S  -  B")  (64) 


where  =  -E  +  iG 

=  complex  conjugate  of 
K2  =  E  -  iH 

K2  =  complex  conjugate  of  K2 


Of  greater  significance  than  the  complex  potential  function  is  its 
derivative  W(S),  which  is  the  complex  velocity  in  the  S-plane. 


dF(S) 


(S  -  A')  (S  -  A")  (S  -  B')  (S  -  B") 


To  get  the  complex  velocity  in  other  planes,  the  chain  rule  is  applied. 
For  the  Z  plane,  W(Z)  is  obtained  as  follows: 


=  dF(Z)  _  dF (S)  dS  dQ 
dZ  dS  dQ  dZ 


(66) 


or  for  the  T-plane: 
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wm  =  dF(T)  =  dF (S)  dS  dQ  dZ_  dZj_ 

K  }  dT  dS  dQ  dZ  dZ'  dT 

Taking  the  required  derivatives  of  the  equations  representing  the 
various  transformations  yields: 


dS  S2  -  d2 
dQ  =  Q2  -  r2d2 


(68) 


2 

dQ  =  Q  -  Q(A  +  B)  +  AB 
dZ  A  -  B 


(69) 


dZ 

dY7  =  a  constant  which  depends  on  which  radial 

cascade  is  desired  (70) 


dZ'  -Z' 
dT  6 


(71) 


The  preceding  equations  complete  all  the  equations  which  are  required 
to  calculate  the  velocity  at  any  point  in  the  various  planes.  It  should 
be  noticed  that  to  calculate  the  derivatives  in  Equations  (68),  (69)  and 
(71),  the  complex  coordinates  of  the  point  under  consideration  are  required 
in  the  S,  Q  and  Z'  planes.  Thus,  to  calculate  the  velocity  at  a  point  in 
the  T-plane,  for  example,  the  corresponding  points  in  the  other  planes  must 
first  be  determined  using  the  transformations  of  Equations  (54),  (55)  and 
(62).  After  the  coordinates  of  these  points  are  known,  the  velocity  at  the 
original  point  in  the  S-plane  is  calculated  using  Equation  (65).  This 
complex  velocity  is  then  multiplied  by  the  derivatives  calculated  from 
Equation  (68)  through  (71)  yielding  the  desired  velocity  at  the  point 


in  the  T-plane. 
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Before  concluding  this  section  on  calculating  an  exact  potential 
flow  solution  for  an  analytically  derived  cascade  geometry,  the  constants 
K^,  K^,  and  remain  to  be  specified.  These  constants  are  chosen  so 
as  to  yield  the  desired  inlet  or  average  flow  condition  and  also  to  satisfy 
the  Kutta  condition.  Since  and  are  complex  conjugate  of  and  K_, 
respectively,  it  is  only  necessary  to  determine  and  K^.  The  desired  flow 
conditions  are  known  or  specified  far  upstream  and  far  downstream  of  the 
cascade  in  the  Z-plane.  Thus,  it  is  desirable  to  express  in  terms  of 
W(z)  as  xz  approaches  plus  infinity  and  in  terms  of  W(Z)  as  approaches 
negative  infinity.  The  points  A  and  B  go  to  plus  infinity  and  minus 
infinity,  respectively,  during  the  transformation  from  the  Q-plane  to 
the  Z-plane.  Thus,  taking  the  limit  of  W(z)  as  S  approaches  A'  and  Q 
approaches  A  yields  the  limit  of  W(Z)  as  x^  approaches  plus  infinity. 

From  this  limit  the  following  equation  results: 

W(Z)  =  K  =  V  -  iV  ,  at  x  =  +  »  .  (72) 

i  x  y  z 

Similarly,  taking  the  limit  of  W(Z)  as  S  approaches  B'  and  Q  approaches 
B  yields: 


W(Z)  =  -  K  =  V  -  iV  ,  at  x  =  -  «  .  (73) 
2  x  y  z  v  ' 

The  velocity  components  in  Equations  (72)  and  (73)  can  be  expressed 
in  terms  of  the  magnitude  of  the  average  cascade  velocity  U,  the  cascade 
average  flow  angle  a  and  the  upwash  and  downwash  velocity  created  by  the 
circulation  on  the  airfoils. 


K,  =  U[cos  a  -  i(sin  a  -  V,  )] 
i  an 


(74) 
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^  =  U[-cos  a  +  i(sin  a  +  V  )]  (75) 

(V  and  V,  are  nondimensionalized  by  U) 
up  dn 

In  order  for  the  velocities  and  coefficients  calculated  with  the  exact 
solution  to  match  the  data  from  the  Douglas  Neumann  program,  and 
are  nondimensionalized  by  U  yielding: 

K1  =  cos  a  -  i(sin  a  -  V  )  (76) 


K_  =  -cos  a  +  i(sin  a  +  V  ) 
2  up 


(77) 


(Note:  The  equations  for  the  exact  solution  in  Reference  [1]  are 


nondimensionalized  by  the  magnitude  of  the  inlet  velocity  V  ) 

!  I 


By  definition,  V  is  equal  to  V,  .  Thus,  if  a  is  specified,  V 

up  n  dn  r  up 


and  Vdn  constitute  one  unknown  which  can  be  solved  for  by  requiring  the 
trailing  edge  of  the  airfoils  to  be  a  stagnation  point.  However,  rather 


than  arbitrarily  selecting  a,  V  and  were  expressed  in  terms  of 


the  lift  coefficient  which  was  given  the  arbitrary  value  of  1.75 


for  the  current  check  cases.  With  C  specified,  a  must  be  determined 

Li 


such  that  the  Kutta  condition  is  satisfied.  The  expressions  for  and 


K2  in  terms  of  a  and  CL  are: 


K1  =  cos  a  -  i (sin  a  -  Tjj p)  , 


(78) 


K2  *  -  cos  a  +  i(sin  a  +  ^gp)  » 


(79) 


1  c n  j  *  -i  .  •  cascade  spacing 

where  SP  =  dimensionless  cascade  spacing  =  - ; - L, - 

chord 


Inserting  the  preceding  expressions  for  and  along  with  their  complex 
conjugates  into  Equation  (65)  provides  the  necessary  equation  for 
determining  a.  The  complex  velocity  W(S)  must  be  zero  at  the  point  on 


i 
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the  cylinder  where  it  intersects  the  xg  axis,  since  this  point  becomes 
the  trailing  edge  of  the  airfoils  in  the  Z-plane.  Solving  for  a  with 
a  simple  interval  halving  technique  yields  a  =  -0.86154  degrees. 

The  determination  of  a  completes  the  explanation  of  the  exact 
analytical  cascade  solution.  Important  parameters  for  the  cascade  in 
the  Z-plane  are  summarized  below: 

SP  =  0.79675 
CL  =  1.75 

a  =  -0.86154  degrees 
cXj  =  28.10800  degrees 
a„  =  -29.43192  degrees 

£i 

Aa  =  57.53992  degrees 

The  values  given  above  for  the  angles  are  also  true  for  the  circular 
cascades  since  angles  are  preserved  by  conformal  transformations. 

4 . 3  Comparisons  of  the  Program  Results  With  Exact  Analytical  Solutions 

The  first  comparison  of  the  program  results  with  an  exact  analytical 
solution  is  made  for  the  linear  cascade  illustrated  in  the  Z-plane  of 
Figure  7.  The  pressure  distributions  on  these  airfoils  obtained  from 
the  cascade  program  and  the  exact  solution  are  shown  in  Figure  9.  It  is 
seen  from  Figure  9  that  the  agreement  in  the  data  is  excellent,  as  it 
should  be  if  the  program  is  functioning  properly.  Near  the  center  portion 
of  the  airfoil,  the  differences  between  the  two  results  are  the  largest. 
For  the  computer  run,  the  points  representing  the  airfoil  were  more 
closely  spaced  near  the  leading  and  trailing  edges  of  the  airfoil  in 
order  to  adequately  describe  the  airfoil  shape.  More  points  in  the 
center  region  of  the  airfoil  should  improve  the  results  in  this  area. 

Near  the  leading  and  trailing  edges  where  a  large  number  of  points  were 
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—  EXACT  analytical  solution 

O  CASCADE-PROGRAM  SOLUTION 


CL  =  1.75 


dj  =  28.108 

aE  =  -29.432 
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Figure  9.  A  Comparison  of  the  Pressure  Distributions  Obtained  From  the 
F.xact  Solution  and  the  Computer  Program  for  the  Analytically 
Derived  Airfoil  Shape  in  a  Linear  Cascade.  The  Cascade 
Geometry  is  Illustrated  in  the  Z-Plane  of  Figure  7 
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used,  the  steep  pressure  gradients  are  very  accurately  calculated  by  the 
cascade  program. 

The  cascade  solution  yielding  the  pressure  distribution  shown  in 
Figure  9  was  obtained  by  representing  the  airfoil  with  99  points.  As  just 
stated,  more  points  should  yield  a  more  accurate  solution,  but  even  with 
99  points  the  results  are  probably  within  required  engineering  accuracy. 
Most  airfoil  shapes  employed  in  engineering  practice  will  not  have  the 
large  amount  of  blade  curvature  present  in  the  analytical  profile.  For 
this  reason,  100  or  less  points  will  usually  adequately  represent  any 
airfoil  of  engineering  significance.  A  means  of  checking  if  a  given 
number  of  points  is  sufficient  to  represent  an  airfoil  is  to  change  the 
number  of  points  and  determine  if  the  solution  is  significantly  altered. 

It  is  recalled  that  the  technique  employed  in  the  Douglas  Neumann  program 
will  yield  results  which  approach  the  exact  solution  as  the  number  of 
points  approaches  infinity. 

In  order  to  verify  that  the  program  functions  properly  for  each  of  the 
various  input  options,  several  cases  were  run  for  each  cascade  with  the 
different  cascade  operating  parameters  given  as  input.  In  addition,  the 
options  of  running  a  case  as  a  subcase  of  the  previous  case  or  directing 
the  program  to  go  directly  to  the  combination  equations  were  checked. 

Using  the  linear  cascade  as  an  example,  CT  was  the  cascade  operating 
parameter  given  as  input  for  the  first  case.  For  the  second  case,  a  was 
the  input  parameter  and  the  case  was  run  as  a  subcase  of  the  first  case. 

For  the  third  and  fourth  cases,  Aa  and  ct^  were  the  input  parameters, 
respectively,  and  the  program  was  directed  to  go  directly  to  the  combina¬ 
tion  equations.  The  pressure  distributions  calculated  by  the  program  for 
each  of  the  four  cases  were  found  to  be  nearly  identical.  Comparisons 
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of  the  cascade  parameters  calculated  by  the  program  for  the  four  cases 
along  with  the  values  obtained  from  the  exact  analytical  solution  are 
presented  in  Table  1. 


TABLE  1 

A  Comparison  of  Cascade  Parameters  Obtained 
From  the  Program  for  the  Various  Input  Parameters 


Case 

Input 

Values  Obtained  From  the  Cascade  Program 

// 

Parameter 

CL 

a 

Aa=araE 

ai 

“E 

1 

CL 

-0.7712 

57.5388 

28.1768 

-29.3620 

2 

a 

1.7451 

57.4040 

28.0392 

-29.3648 

3 

Aa 

1.7503 

-0.7645 

28.1871 

-29.3618 

4 

aI 

1.7475 

-0.8165 

57.4714 

-29.3634 

J 

Values  From  the 
Exact  Solution 

1.75 

-0.8615 

57.5399 

28.1080 

-29.4319 

From  the  fact  that  the  pressure  distributions  for  the  four  cases  were 
nearly  identical  and  from  the  results  shown  in  Table  1,  it  is  apparent 
that  the  various  input  and  running  options  in  the  program  function  properly 
The  linear  cascade  shown  in  the  Z-plane  of  Figure  7  will  also  be 
used  for  an  example  of  the  verification  of  the  program  results  for  off-body 
points.  The  locations  of  the  off-body  points  are  marked  with  "x's"  in 
Figure  7.  Table  2  lists  the  coordinates  of  the  off-body  points  along 
with  the  velocity  components  at  these  points  obtained  from  both  the 
cascade  program  and  the  exact  analytical  solution.  The  values  shown 
in  the  table  verify  that  the  program  properly  calculates  the  velocity 
components  at  off-body  points. 
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TABLE  2 

A  Comparison  of  Velocity  Components  at  Off-Body  Points 
Obtained  From  the  Cascade  Program  and  the  Exact  Analytical  Solution 


Off-Body  Point 
Coordinates 

Cascade  Program 
Results 

Exact  Analytical 
Solution 

x/ c 

y/c 

V  /u 

X 

V  /U 

y 

v  /u 

X 

V  /u 
y 

-0.7424 

0.3984 

0.9663 

0.5518 

0.9683 

0.5502 

-0.7424 

0.0 

1.0239 

0.5195 

1.0252 

0.5180 

-0.7424 

-0.3984 

0.9663 

0.5518 

0.9683 

0.5502 

-0.4924 

0.1992 

0.9756 

0.3831 

0.9762 

0.3829 

-0.4924 

0.0 

1.1106 

0.4345 

1.1110 

0.4333 

-0.4924 

-0.1992 

1.2289 

0.6373 

1.2285 

0.6344 

0.0 

0.3011 

0.9510 

0.0088 

0.9526 

0.0086 

0.0 

0.1188 

1.0827 

-0.0028 

1.0830 

-0.0031 

0.0 

-0.0635 

1.2305 

-0.0133 

1.2306 

-0.0138 

0.5077 

0.1992 

0.9427 

-0.4307 

0.9404 

-0.4314 

0.5077 

0.0 

1.0824 

-0.4592 

1.0809 

-0.4601 

0.5077 

-0.1992 

1.1879 

-0.6175 

1.1868 

-0.6188 

0.7577 

0.3984 

0.9806 

-0.5822 

0.9790 

-0.5840 

0.7577 

0.0 

1.0190 

-0.5445 

1.0174 

-0.5461 

0.7577 

-0.3984 

0.9806 

-0.5822 

0.9790 

-0.5840 

In  Figure  10,  ti  e  pressure  distributions  are  presented  for  the 
airfoils  in  the  circular  cascade  with  the  radial  flow  inward.  This 
cascade  is  illustrated  in  the  T-plane  of  Figure  7.  As  seen  from  Figure  10, 
the  results  from  the  exact  analytical  solution  and  the  cascade  program  are 
nearly  identical.  Again,  the  differences  in  the  results  mould  be  reduced 
by  increasing  the  number  of  points. 

It  is  interesting  to  compare  the  shape  of  the  pressure  distribution 
for  the  airfoils  in  the  linear  cascade  with  the  pressure  distribution 
for  the  airfoils  in  the  circular  cascade.  From  Figure  9,  it  is  seen 
that  the  pressure  distribution  for  the  linear  cascade  has  two  nearly 


I 
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Figure  10.  A  Comparison  of  the  Pressure  Distributions  Obtained  From  the 
Exact  Solution  and  the  Computer  Program  for  the  Circular 
Cascade  Illustrated  in  the  T-Plane  of  Figure  7 
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symmetrical  peaks.  The  pressure  distribution  shown  in  Figure  10  for 
the  circular  cascade  with  the  radial  flow  inward  has  one  large  peak 
and  one  small  peak.  The  airfoil  geometries  in  the  two  cascades  are 
a  little  different  due  to  the  transformation  from  the  Z-plane  to  the 
T-plane.  However,  the  major  effect  causing  the  pressure  distributions 
to  have  different  shapes  is  the  change  in  radius  for  the  circular  cascade 
The  trailing  edge  of  the  airfoil  in  the  circular  cascade  is  at  a  smaller 
radius  than  the  leading  edge.  For  this  reason,  the  velocity  near  the 
trailing  edge  will  be  greater  than  the  velocity  near  the  leading  edge, 
resulting  in  the  trailing  edge  peak  to  be  greater  than  the  peak  in  the 
pressure  distribution  near  the  leading  edge. 

Figure  11  shows  the  comparison  of  the  pressure  distributions  for 

the  last  type  of  cascade,  the  circular  cascade  with  the  radial  flow 

outward.  This  cascade  is  illustrated  in  the  T-plane  of  Figure  8.  In 

Figure  8,  notice  the  orientation  of  the  airfoil  in  the  local  coordinate 

system  used  for  the  program  input.  Since  the  airfoil  is  inverted  in  the 

local  coordinate  system,  the  airfoil  is  also  shown  this  way  in  Figure  11. 

In  addition,  the  scale  for  C  is  inverted  from  the  previous  scales  to 

P 

match  the  airfoil.  For  the  circular  cascade  with  the  radial  flow  outward 
it  is  again  seen  that  the  pressure  distributions  computed  with  the  exact 
analytical  solution  and  the  rewritten  version  of  the  Douglas  Neumann 
program  are  nearly  identical. 

The  pressure  distribution  shown  in  Figure  11  is  similar  to  the 
pressure  distribution  in  Figure  10,  in  that  it  has  one  large  peak  and 
one  small  peak.  However,  for  the  circular  cascade  with  the  radial  flow 
outward,  the  large  peak  is  near  the  leading  edge,  since  the  leading 
edge  is  at  the  smaller  radius  where  the  velocity  is  Larger. 
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igure  11.  A  Comparison  of  the  Pressure  Distributions  Obtained  From  the 
Kxact  Solution  and  the  Computer  Program  for  the  Circular 
Cascade  Illustrated  in  the  T-Plane  of  Figure  8 
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For  both  types  of  circular  cascades,  the  various  input  options, 
program  running  options,  and  the  off-body  point  results  were  all 
verified  in  the  same  manner  as  explained  for  the  linear  cascade.  The 
three  types  of  single  cascades  were  then  modeled  as  multiple  cascades  to 
verify  the  program  results  for  multiple  cascade  analyses.  The  pressure 
distributions  for  the  different  cascades  of  a  multiple  cascade  model 
were  found  to  be  nearly  identical  to  each  other  and  to  the  results 
for  the  single  cascade.  The  flow  angles  and  off-body  point  results 
were  also  found  to  be  correct  for  multiple  cascade  analyses. 

4 . 4  Verification  of  Results  for  Multiple  Cascades  of  Lifting  and 

Nonlifting  Bodies 

A  special  case  of  a  multiple  cascade  analysis  results  when  some 
cascades  have  lifting  bodies  and  others  have  nonlifting  bodies.  To 
verify  that  the  current  version  of  the  Douglas  Neumann  program  could 
handle  this  special  case,  results  from  the  current  version  of  the  program 
were  compared  to  results  published  for  the  original  program  in 
Reference  f 1 ] .  The  cascades  analyzed  and  the  resulting  pressure  distri¬ 
butions  are  shown  in  Figure  12.  The  first  and  third  cascades  consist  of 
cylinders  with  no  circulation.  The  second  cascade  consists  of  cylinders 
which  have  circulation.  In  the  computer  program,  the  circulation 
resulted  from  specifying  the  rear  stagnation  point  (trailing  edge) 
to  be  at  -30°  on  the  cylinder.  As  seen  from  Figure  12,  the  pressure 
distributions  from  the  two  programs  are  nearly  identical,  which  verifies 
the  current  version  of  the  program  can  correctly  analyze  the  special 
case  of  multiple  cascades  with  lifting  and  nonlifting  bodies.  (Note: 

The  first  and  third  cascades  do  have  lift  forces  on  them  due  to  induced 
effects  of  the  second  cascade.  However,  they  do  not  produce  lift  on 


their  own.) 


Figure  12.  A  Comparison  of  Pressure  Distributions  for  Cascades  of  Lifting 

and  Non-Lifting  Bodies  Obtained  From  the  Original  and  the  Current 
Version  of  the  Douglas  Neumann  Cascade  Program 
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APPENDIX  A 

A  Users  Guide  for  Running  the  Program 
A. 1  Explanation  of  Program  Input 

The  main  body  of  this  report  has  given  the  theory  behind  the  program, 
has  defined  the  input  and  output  variables,  and  has  discussed  how  various 
types  of  cascades  can  be  modeled.  A  basic  understanding  of  the  preceding 
material  should  enable  the  program  to  be  used  and  the  results  understood 
with  little  difficulty.  This  appendix  will  give  detailed  instructions 
on  how  to  input  the  required  data  to  the  program.  In  the  latter  sections 
of  this  appendix,  sample  input  and  sample  output  are  presented. 

Before  describing  the  program  input,  attention  is  called  to  Figures  A1 
ad  A2  which  illustrate  the  geometric  parameters  involved  in  the  input. 
Figure  A1  illustrates  two  linear  cascades  and  Figure  A2  illustrates 
two  circular  cascades.  Input  for  one  cascade  or  more  than  two  cascades 
follow  the  same  format.  With  regard  to  circular  cascades,  it  should  be 
pointed  out  that  except  for  two  parameters,  the  input  for  circular 
cascades  with  the  radial  flow  inward  and  the  input  for  circular  cascades 
with  the  radial  flow  outward  are  identical.  Of  the  two  parameters  which 
differ,  one  tells  the  program  the  flow  direction  and  the  other  specifies 
the  location  of  the  trailing  edge  when  the  radial  flow  direction  is 
outward, 

For  each  input  card  required  by  the  program,  the  title  of  the 
card  will  be  followed  by  the  variables  which  appear  on  the  card  and  their 
required  format.  After  the  list  of  variables,  a  description  of  each 
variable  will  be  given.  The  following  sequence  of  cards  can  be  repeated 
several  times  in  order  to  run  several  cases  at  once.  The  program  will 
continue  to  process  runs  until  a  code  on  the  second  card  tells  the 
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program  to  terminate.  The  various  options  available  for  analyzing  a 
special  case  of  the  previous  run  will  be  explained  as  the  appropriate 
input  variables  are  encountered. 

Heading  Card:  HEDR  FORMAT  (20A4) 

HEDR  -  A  heading  which  describes  the  run  and  is  read  as  alphanumeric 
characters.  This  heading  is  printed  at  the  top  of  each  page 
of  output. 

Control  Card:  NB,  FLG02,  FLG03,  FLG04,  FLG05,  FLG06,  FLG07 ,  FLG08, 

FLG09,  FLG10,  FLG11 ,  FLG12  FORMAT  (1215) 

NB  -  number  of  bodies  or  cascades 

FLG02  -  A  nonzero  integer  if  the  flow  is  to  be  determined  at  points 
off  of  the  bodies. 

As  stated  in  the  main  body  of  this  report,  a,  Aa,  a^.,  or  C^  can  be 
used  in  the  combination  equations  for  specifying  the  operating  condition 
of  the  cascades.  FLG03  through  FLG06  specify  which  input  parameter  is 
used . 

FLG03  -  A  nonzero  integer  if  a  is  to  be  used  in  the  combination  equations. 

FLG04  -  A  nonzero  integer  if  Aa  is  to  be  used  in  the  combination  equations 

FLG05  -  A  nonzero  integer  if  a^.  is  to  be  used  in  the  combination  equations 

FLG06  -  A  nonzero  integer  if  C  is  to  be  used  in  the  combination  equations 

1j 

C^  cannot  be  used  as  input  for  circular  cascades. 

FLG07  -  A  nonzero  integer  if  the  matrices  of  influence  coefficients  (A  and 
B)  and  the  onset  flows  are  to  be  printed.  The  A  matrix  represents 
the  velocity  components  induced  normal  to  the  surface  and  the 
B  matrix  represents  the  velocity  components  induced  tangent  to 
the  surface.  A  and  B  are  large  matrices  and  normally  are  not 
printed  (FLG07  =  0). 


Schematic  of  Two  Linear  Cascades  Illustrating  the  Geometric  Input 
Parameters  for  the  Program 
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Figure  A. 2.  Schematic  of  Two  Circular  Cascades  Illustrating  the  Geometric 
Input  Parameters  for  the  Program.  The  Cascade  Geometric 
Input  are  Independent  of  the  Radial  Flow  Direction 
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FLG08 


FLG09 


FLG10 


FLG11 


FLG12 


A  nonzero  integer  if  the  program  is  to  go  directly  to  the  combination 
equations  using  the  basic  solutions  of  the  previous  case.  This 
option  can  be  used  only  if  the  cascade  geometry  is  identical  to 
the  previous  case  (i.e.  only  the  flow  angle  or  CT  is  changed). 
Off-body  point  data  input  is  independent  of  FLG08. 

A  nonzero  integer  if  the  cascades  under  consideration  are  circular 
cascades . 

For  circular  cascades,  FLG10  is  a  nonzero  integer  if  the  radial 
flow  is  outward  and  FLG10  is  zero  if  the  radial  flow  is  inward. 

A  nonzero  integer  if  the  first  cascade  geometry  is  to  be  repeated 
for  a  multiple  cascade  analysis.  This  option  was  included  as  an 
aid  in  inputting  data  when  modeling  the  cascades  requires  one 
geometry  to  be  repeated  with  equal  spacing  between  each  body 
(as  with  the  wicket  gates  in  the  example  illustrated  in  Figure  6) . 

A  nonzero  integer  if  the  program  is  to  be  terminated.  The 
program  will  continue  to  try  to  process  additional  runs  until 
terminated.  Thus,  terminating  the  program  at  the  end  of  a  run 
requires  a  heading  card,  which  can  be  left  blank,  and  a  control 
card  with  FLG12  a  nonzero  integer. 


i 

i 


J 
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Cascade  Parameter  Card:  SP,  CL,  ALPHA,  FALPHA,  DALPHA,  SUBKS  FORMAT  (5F10.0,  15) 
SP  -  Cascade  Spacing 

For  linear  cascades,  SP  is  the  dimensionless  cascade  spacing  nondimen- 
sionalized  by  the  chord  of  one  of  the  cascades.  The  chord  used  to 
nondimensionalize  SP  must  also  be  used  to  nondimens ionalize  the  chords 
of  all  the  cascades  in  order  to  keep  the  relative  size  and  spacing 
correct.  From  this  point  on,  the  chord  of  the  body  used  to  nondimen¬ 
sionalize  SP  will  be  referred  to  as  the  reference  chord  c 

ref 

For  circular  cascades,  SP  is  the  angular  spacing  in  degrees. 

Only  one  of  the  parameters  CL,  ALPHA,  FALPHA,  DALPHA  needs  to  be 
specified,  depending  on  which  one  is  to  be  used  in  the  combination 
equations.  If  not  used  in  the  combination  equations,  any  value  read  is 
ignored . 

CL  -  Lift  coefficient.  CL  can  only  be  input  for  linear  cascades. 

ALPHA  -  Average  flow  angle  a. 

FALPHA  -  Inlet  flow  angle  a^.. 

DALPHA  -  Change  in  flow  angle  through  the  cascades,  Aa. 

SUBKS  -  A  nonzero  integer  if  the  unmodified  coordinates  of  all  the 
bodies  of  the  previous  case  are  to  be  used  for  the  present 
case.  The  chords,  stagger  angles,  position  of  the  cascades, 
etc.  may  change.  For  a  subcase,  the  program  reads  the  coordinate 
data  from  files  generated  in  the  previous  run,  but  then  proceeds 
just  as  if  the  coordinate  data  had  been  read  from  cards. 

If  FLG08  is  nonzero  (i.e.  the  program  is  going  directly  to  the  combination 
equations),  the  cascade  parameter  card  completes  the  required  input  describing 
the  cascades.  The  next  input  will  be  the  off-body  point  data,  if  the  flow 
at  off-body  points  is  to  be  determined  for  this  case. 
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Similar  Cascade  Card:  NCSG  FORMAT  (15) 

NCSG  -  The  number  of  cascades  with  geometry  similar  to  cascade  1, 
including  cascade  1. 

This  card  is  only  required  if  several  cascades  have  the  same  geometry  and 
the  coordinates  of  the  similar  cascades  are  to  be  generated  internally 
subsequent  to  the  input  of  the  coordinates  for  the  first  cascade  (i.e. 
this  card  is  not  required  unless  FLG11  is  nonzero). 

Cascade  Geometry  Card:  CHORD,  CCANG,  XP,  YP,  RPADDX ,  ADDY,  NP,  NLF, 

NLE  FORMAT  (6F10.0,  315) 

CHORD  -  Dimensionless  chord  of  the  body.  CHORD  should  be  nondimensionalized 
by  the  chord  of  one  of  the  cascades  c  For  a  single  cascade 

analysis,  CHORD  =  1. 

CCANG  -  Stagger  angle  £  of  the  airfoils  in  the  cascade,  as  illustrated 
in  Figures  A1  and  A2. 

XP,  YP  -  Local  coordinates  of  a  point  in  or  on  the  body  which  is  used  to 
locate  the  body  in  the  global  coordinate  system.  The  body  is 
also  rotated  around  the  point  (XP,  YP)  to  obtain  the  desired 
stagger  angle.  In  addition,  the  force  and  moment  coefficients 
are  calculated  at  this  point.  XP  and  YP  must  be  nondimensionalized 
by  the  chord  of  the  body  they  are  for. 

RPADDX  -  For  circular  cascades,  RPADDX  is  the  radius  in  the  global  coordinate 
system  at  which  the  point  (XP,  YP)  is  located,  as  illustrated  in 
Figure  A2. 

For  linear  cascades,  RPADDX  is  the  x  location  in  the  global  coordinate 
system  at  which  the  point  (XP,  YP)  is  located,  as  illustrated  in 
Figure  Al. 

RPADDX  must  be  nondimensionalized  by  c 


— 
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ADDY  -  For  circular  cascades,  ADDY  is  the  angular  location  in  the  global 

coordinate  system  of  the  point  (XP,YP)  for  one  of  the  bodies  in  the 
cascade.  ADDY  is  illustrated  in  Figure  A2  and  must  be  expressed  in 
degrees . 

For  linear  cascades,  ADDY  is  the  y  location  in  the  global  coordinate 

system  of  the  point  (XP,YP)  for  one  of  the  bodies  in  the  cascade. 

ADDY  is  illustrated  in  Figure  A1  and  must  be  nondimensionalized 

by  c  ... 
ref 

NP  -  Number  of  points  which  will  be  used  to  represent  the  body. 

NLF  -  A  nonzero  integer  if  the  cascade  is  nonlifting.  If  lifting 
and  nonlifting  bodies  are  both  present  in  a  multiple  cascade 
analysis,  the  data  for  the  lifting  bodies  must  be  read  in  first. 

NLE  -  For  circular  cascades  with  the  radial  flow  outward,  NLE  is  the 
number  of  the  point  which  will  be  the  true  trailing  edge  of  the 
airfoil.  NLE  is  not  required  for  the  other  types  of  cascades, 
because  the  location  of  the  trailing  edge  is  known  from  the 
starting  point  of  the  data. 

Body  Coordinate  Cards:  X,  Y  FORMAT  (2F10.0) 

X,  Y  -  Coordinates  of  points  on  the  body  in  the  local  coordinate  system 
x'-y'  which  describe  the  shape  of  the  body,  x  and  y  must  be 
nondimensionalized  by  the  chord  of  the  body.  There  are  NP 
body  coordinate  cards  and  the  data  must  start  at  the  points 
labeled  in  Figures  A1  and  A2.  For  linear  cascades  and  circular 
cascades  with  the  radial  flow  inward,  the  starting  point  is  the 
trailing  edge  of  the  airfoil.  The  first  point  must  be  repeated 


as  the  last  point. 
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Unless  the  case  being  run  is  a  subcase  or  similar  cascade  geometries  are 
being  generated  internally,  the  cascade  geometry  card  and  the  body  coordinate 
cards  are  repeated  for  each  cascade  of  a  multiple  cascade  analysis.  If  the 
case  being  run  is  a  subcase  (SU3KS  is  nonzero),  only  the  cascade  geometry 
cards  are  required  and  the  body  coordinate  data  will  be  read  from  files. 

If  similar  cascades  are  being  generated  internally  (FLG11  is  nonzero) , 
the  cascade  geometry  card  and  the  body  coordinate  cards  are  required  for 
the  first  cascade  and  for  each  different  cascade  after  the  first  NCSG 
cascades.  For  the  cascades  generated  internally,  the  parameters  CHORD, 

CCANG,  XY,  YP,  RPADDX,  NP,  NLF  and  NLE  are  all  set  equal  to  the  values 
for  the  first  cascade.  ADDY  is  indexed  by  the  value  SP/NCSG  for  each 
cascade  generated  internally. 

Off-Body  Point  Data 

Coordinates  of  off-body  points  where  the  flow  is  to  be  analyzed 
can  be  read  in  or  generated  internally  after  providing  initial  information 
to  the  program.  The  off-body  points  can  be  grouped  in  sets  in  order  to 
pro’-ide  flexibility  in  the  generation  and  distribution  of  the  points. 

None  of  the  following  cards  are  required  if  FLG02  is  zero. 

Number  of  Sets  Card:  NSETS  FORMAT  (15) 

NSETS  -  Number  of  sets  of  off-body  points  to  be  read  in  or  generated 
internally. 

Set  Initial  Data  Card:  NCR,  NFS,  XR,  YT1,  DELYT  FORMAT  (215,  3F10.0) 

NCR  -  A  nonzero  integer  if  the  off-body  points  of  this  set  are  to  be 
read  from  cards 


-  Zero  if  the  off-bodv  points  of  this  set  are  to  be  generated  internal] y. 
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NPS  -  Number  of  points  in  the  set 

If  NCR  is  a  nonzero  integer,  XR,  YT1  and  DELYT  do  not  need  to  be  specified. 

XR  -  For  linear  cascades,  XR  is  the  x  coordinate  of  the  off-body  points 
to  be  generated  in  this  set. 

For  circular  cascades,  XR  is  the  r  coordinate  of  the  off-body 

points  to  be  generated  in  this  set. 

XR  must  be  nondimensionalized  by  c  .. 

ref 

YT1  -  For  linear  cascades,  YT1  is  the  y  coordinate  of  the  first  point 

in  the  set.  YT1  must  be  nondimensionalized  by  c  _. 

ref 

For  circular  cascades,  YT1  is  the  0  coordinate  of  the  first  point 
in  the  set,  expressed  in  degrees. 

DELYT  -  For  linear  cascades,  DELYT  is  the  y  increment  by  which  the 

remainder  of  the  points  in  the  set  will  be  generated.  DELYT 

must  be  nondimensionalized  by  c  , . 

ref 

For  circular  cascades,  DELYT  is  the  0  increment  by  which  the 
remainder  of  the  points  in  the  set  will  be  generated.  DELTY 
must  be  expressed  in  degrees. 

Off-Body  Coordinate  Cards:  XOB,  YOB  FORMAT  (2F10.0) 

If  the  off-body  points  are  to  be  read  in  for  this  set  (NCR  is  nonzero), 

NPS  off-body  point  coordinate  cards  are  required. 

XOB,  YOB  -  For  linear  cascades,  XOB,  YOB  are  the  x,  y  coordinates  of  the 
off-body  points  nondimensionalized  by  c 

ref 

For  circular  carcades,  XOB,  YOB  are  the  r,  0  coordinates, 
respectively,  of  the  off-body  points.  XOB  must  be  nondimensionalized 
by  c  ^  and  YOB  must  be  expressed  in  degrees. 
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The  set  initial  data  card  and  the  off-body  point  coordinate  cards  when 
appropriate  are  required  for  each  set. 

The  off-body  point  data  completes  the  description  of  the  various  types 
of  input  for  the  program.  Data  for  the  next  case  can  follow  the  present 
case  starting  with  a  new  heading  card.  As  previously  stated,  the  program 
can  be  terminated  by  supplying  a  blank  heading  card  and  making  FLG12  a 
nonzero  integer  on  the  control  card. 


I 
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APPENDIX  A. 2 

Sample  Program  Input: 

Analytic  Cascade  Profile  C — Multiple  Rectangular  Cascades 

The  following  sample  input  data  represent  the  linear  cascade  and 
off-body  points  shown  in  the  Z-plane  of  Figure  7.  The  cascade  is  modeled 
as  two  cascades  to  illustrate  a  multiple  cascade  analysis.  The  coordinates 
for  the  second  cascade  are  generated  internally  (FLG11  =  nonzero  integer 
and  NCSG  =  2),  thus,  only  one  set  of  coordinate  data  points  are  given. 
Relative  to  the  airfoils,  the  off-body  points  are  the  same  as  those  shown 
in  Figure  7  and  Table  2.  However,  since  the  leading  edge  of  the  reference 
airfoil  was  positioned  at  the  origin  of  the  global  system,  the  off-body 
points  reflect  this  position  and  are  translated  from  those  shown  in  the 
figure  and  table. 

Data  are  given  for  four  consecutive  runs  which  demonstrate  the  programs 

various  input  options.  The  first  run  has  C  as  the  given  cascade  operating 

i-* 

parameter.  (Note:  C  =  2(1.75)  =  3.5,  since  the  cascade  is  modeled  as 

Li 

two  similar  cascades.)  Following  the  data  for  the  first  run,  the  data  for 
the  second  run  directs  the  program  to  go  directly  to  the  combination 
equations  using  Aa  as  the  cascade  operating  parameter.  The  third  run 
is  a  subcase  of  the  previous  run  (coordinate  data  are  read  from  a  file) 
and  01^  is  the  given  cascade  parameter.  The  final  run  again  directs  the 
program  to  go  directly  to  the  combination  equations,  this  time  with  a 
as  the  cascade  operating  parameter.  The  data  for  these  four  run  options 
will  enable  the  user  to  verify  that  the  program  is  operating  properly. 

The  output  for  the  first  run  only  is  given  in  Section  A. 3.  The  output 
for  the  subsequent  runs  should  be  similar. 
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APPENDIX  A. 3 


Sample  Program  Output  Using  Data  in  Appendix  A. 2. 

Only  the  Output  Which  Uses  CT  in  the  Combination 
Equations  is  included. 


INPUT  PARAMETERS 


A.IAirriC  CASCACE  PROFILE  C — ftfJLTIPLR  HEC1A NGOIAF  CASCADES 
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APPENDIX  B 

B. 1  Data  Storage  and  Equation  Solution  Procedures 

Appendix  B  is  included  in  this  report  because  the  data  storage  space 
required  by  the  program  is  one  of  the  primary  concerns  when  adapting  the 
program  to  a  particular  computer  system  and/or  a  particular  problem.  A 
discussion  of  the  available  methods  for  solving  the  system  of  equations 
is  included,  because  the  method  used  depends  on  the  size  of  the  problem 
and  the  amount  of  core  storage  available.  Although  the  program  is  quite 
general,  the  purpose  of  this  appendix  is  to  give  the  program  user  an 
understanding  of  the  storage  requirements  and  equation  solution  techniques 
so  that  he  can  modify  the  program  to  meet  his  needs  in  the  most  economical 
and  time  efficient  manner.  For  example,  if  the  user  has  large  problems 
but  a  computer  with  a  small  core;  the  array  dimensions  will  have  to  be  changed 
to  match  the  available  space  and  the  solution  procedure  selected  which  solves 
the  equations  with  the  terms  in  the  equations  read  from  files  or  auxiliary 
storage.  On  the  other  hand,  if  relatively  small  problems  are  run,  the 
entire  matrix  representing  the  system  of  equations  may  be  stored  in  the 
core  and  reduced  in  a  conventional  manner  for  the  solution. 

The  major  portion  of  the  calculations  performed  in  the  program  are 
done  to  solve  the  system  of  simultaneous  equations  which  yield  the  source 
distributions  on  the  surface  of  the  airfoils.  Equation  (11)  in  the  main 
body  of  this  report  represents  the  system  of  equations  in  matrix  form. 

To  aid  in  the  discussion,  Equation  (11)  is  rewritten  below: 
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A11  A12 
A21  A22 


'Su  ^2 


or 


> 

*1 

"l 

'  A2N 

°  2 

- 

(B.l) 


[A]  [a]  =  [0] 


where  [A]  =  Influence  coefficient  matrix 

[a]  =  Source  distribution  matrix  (solution  matrix) 

[0]  =  Onset  flow  matrix 

As  discussed  in  the  main  body  of  the  report,  several  basic  solutions  are 
required  for  each  problem.  Thus,  Equation  (B.l)  must  actually  be  solved 
several  times.  During  each  step  of  the  solution,  the  solution  matrices 
for  each  basic  solution  are  modified;  thus,  the  several  required  basic 
solutions  are  essentially  obtained  simultaneously. 

Before  discussing  the  solution  techniques,  an  important  difference 
should  be  noted  between  the  A  matrix  of  Equation  (B.l)  and  the  matrices  of 
typical  finite  element  and  finite  difference  methods.  The  sources  applied 
to  the  surface  elements  of  the  airfoils  in  the  Douglas  Neumann  method  have 
an  induced  effect  on  every  element.  The  coefficients  in  the  A  matrix  in 
Equation  (B.l)  represent  this  influence;  and  if  there  are  N  elements,  there 
are  N  rows  in  the  matrix  each  with  N  terms.  Typical  finite  difference  and 
finite  element  methods  yield  matrices  with  mostly  zeros,  except  near  the 
diagonal.  For  the  finite  difference  and  finite  element  methods,  only 


-104- 


26  June  1981 
AMY:  cag 


the  nonzero  terms  near  the  diagonal  need  to  be  stored  and  manipulated 

for  a  solution.  In  comparison,  the  terms  in  the  A  matrix  are  all  nonzero, 

2 

and  therefore  N  surface  elements  result  in  N  terms  which  must  be  stored 

2 

and  manipulated  for  the  solution.  With  N  terms,  it  is  easy  to  see  that  for 
a  large  number  of  surface  elements  representing  the  cascades,  storage 
problems  and  long  computer  solution  times  may  be  encountered. 

There  are  currently  three  methods  available  to  solve  the  system  of 
equations  in  the  rewritten  version  of  the  Douglas  Neumann  program.  Two 
of  the  methods  were  in  the  original  Douglas  Neumann  cascade  program  and  a 
third  method  was  developed  by  ARL  which  can  handle  very  large  problems 
and  requires  less  computer  time.  In  the  following  paragraphs,  the  three 
methods  will  each  be  explained  and  their  core  storage  requirements  discussed. 
After  the  three  methods  are  described,  an  example  illustrating  the  computa¬ 
tion  time  required  by  the  three  methods  will  be  given.  The  second  part 
of  the  appendix  then  specifically  discusses  the  array  dimensions  in  the 
program  and  how  these  dimensions  and  the  equation  solution  procedure  can 
be  changed. 

The  first  method  to  be  discussed  for  solving  the  system  of 

simultaneous  equations  was  in  the  original  Douglas  Neumann  cascade 

program.  It  is  a  matrix  reduction  technique  which  requires  that  the 

entire  A  matrix  of  Equation  (B.l)  be  stored  in  an  array  in  the  computer 

core.  For  each  basic  solution  required,  there  is  an  onset  flow  matrix 

[0]  which  must  also  be  in  the  core  if  all  the  basic  solutions  are  to 

be  obtained  with  a  single  reduction  of  the  A  matrix.  If  there  are  N 

elements  and  Ng  basic  solutions  needed,  the  total  array  area  required 

2  9 

by  this  reduction  technique  is  N  +  N  N.  With  an  array  area  of  +  N  N 

s  s 


-105- 


26  June  1981 
AMY : cag 


required,  the  reduction  technique  is  usually  not  practical  for  problems 
with  N  greater  than  two  hundred.  One  hundred  elements  is  usually  sufficient 
to  represent  one  airfoil,  thus  two  cascades  can  easily  be  represented  by 
less  than  two  hundred  elements  and  the  equations  solved  using  the  matrix 
reduction  technique. 

The  idea  behind  the  typical  matrix  reduction  type  of  simultaneous 
equation  solution  is  to  operate  on  the  coefficient  matrix  and  the  matrix 
on  the  right  hand  side  of  the  equation  until  the  coefficient  matrix 
consists  of  ones  on  the  diagonal  and  zeroes  everywhere  else.  When  the 

coefficient  matrix  is  reduced  to  this  form,  the  matrix  on  the  right 

hand  side  of  the  equation  is  the  solution.  Since  each  row  of  the  two 
matrices  represents  an  equation,  rows  can  be  multiplied  by  constants  and 
added  or  subtracted  from  each  other  just  as  equations  are  manipulated. 

As  the  coefficient  matrix  is  reduced,  the  same  operations  are  performed 
on  the  matrix  on  the  right  hand  side  of  the  equation.  However,  the  values 

in  the  matrix  on  the  right  hand  side  do  not  effect  the  reduction  of  the 

coefficient  matrix.  For  this  reason,  in  the  Douglas  Neumann  program  all 
of  the  0  matrices  can  be  operated  on  during  the  reduction  of  the  A  matrix, 
and  all  the  basic  solutions  are  obtained  simultaneously. 

In  a  matrix  reduction  type  of  solution  procedure,  it  is  known  to  be 
good  practice  to  have  one  of  the  larger  values  of  the  row  on  the  diagonal. 
Since  the  rows  are  first  divided  by  the  value  on  the  diagonal,  a  zero  cannot 
be  on  the  diagonal  and  small  values  on  the  diagonal  may  result  in  roundoff 
errors.  In  the  Douglas  Neumann  program,  rather  than  first  rearranging  the 
matrix  for  a  strong  diagonal,  and  equivalent  nrocedure  is  used  where  the 
largest  value  of  the  row  is  used  to  divide  the  other  values  of  row  and  then 


eliminate  other  terms  in  the  same  column.  The  result  is  a  one  in  each  column 
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with  the  remaining  values  zero.  However,  since  the  largest  values  are  not 
necessarily  on  the  diagonal,  the  solution  when  the  reduction  is  completed 
is  not  in  the  correct  order.  The  second  part  of  the  solution  technique 
thus  consists  of  rearranging  the  solution  matrix  so  the  values  are  in  the 
correct  order.  One  improvement  the  author  of  this  report  feels  could  be 
made  would  be  to  use  a  more  efficient  method  of  putting  the  solution  in 
its  proper  order. 

The  second  method  available  for  solving  the  set  of  simultaneous  equations 
was  also  in  the  original  version  of  the  Douglas  Neumann  cascade  program  and 
is  found  in  one  version  of  Subroutine  PART3.  A  method  was  required  which 
could  handle  larger  problems  where  the  A  matrix  was  too  large  to  be  stored 
in  the  core  of  the  computer.  In  the  original  Douglas  Neumann  program  an 
iterative  technique  was  used  which  only  requires  one  row  of  the  A  matrix  to  be 
in  the  computer  core  at  one  time.  A  row  of  the  A  matrix  is  read  from 
auxiliary  storage  and  multiplied  by  the  solution  matrix  [a]  as  it  currently 
appears  ([a]  is  initially  zero).  If  the  solution  matrix  is  correct,  the 
value  resulting  from  the  multiplication  will  equal  the  corresponding  value 
in  the  onset  flow  matrix.  The  difference  in  the  two  values  divided  by 
Ajj  indicates  the  amount  should  be  changed  for  the  current  iteration. 

The  value  of  a ^  is  then  updated  and  while  a  particular  row  of  the  A  matrix 
is  in  the  core  the  process  is  repeated  for  the  other  basic  solutions. 

When  the  calculations  are  complete  for  all  of  the  basic  solutions,  the 
next  row  of  the  A  matrix  is  read  from  auxiliary  storage  and  the  procedure 
is  followed  until  all  the  rows  have  been  treated.  At  the  end  of  a  complete 
iteration  the  maximum  change  in  any  value  of  a  is  compared  to  a  convergence 
criterion  and  the  iterations  continue  until  the  criterion  is  met.  Since 
different  basic  solutions  may  converge  with  a  different  number  of  iterations, 
the  iterations  only  continue  for  those  basic  solutions  which  have  not 
converged. 
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The  main  disadvantage  of  the  iterative  approach  is  that  the  A  matrix 
must  be  read  from  auxiliary  storage  for  every  iteration.  Since  reading 
files  is  a  much  slower  process  than  using  values  in  the  computer  core, 
the  computation  time  for  the  iterative  method  can  be  large. 

For  the  purpose  of  reducing  computation  time  for  large  problems,  a 
third  method  for  solving  the  set  of  simultaneous  equations  was  developed 
at  ARL  during  the  current  rewriting  of  the  Douglas  Neumann  program.  The 
third  method  employs  a  simple  matrix  reduction  technique;  but  for  large 
problems  where  all  the  equations  cannot  fit  in  the  core  simultaneously, 
it  combines  the  use  of  auxiliary  storage  and  computer  core  and  operates 
on  only  a  partial  block  of  the  equations  in  the  core  at  a  given  time. 

The  technique  is  designed  to  maximize  the  use  of  whatever  core  space 
is  available  and,  thus,  the  dimension  of  the  array  used  in  the  reduction 
can  be  made  whatever  size  is  best  for  the  computer  system. 

In  the  third  method,  given  a  work  space  of  a  certain  size,  the  sub¬ 
routine  first  sets  aside  a  portion  of  the  array  sufficient  to  hold  one 
equation.  This  space  is  needed  for  reading  the  writing  data  on  the  files. 
With  the  remaining  space,  the  subroutine  determines  how  many  equations 
can  fit  into  the  array  and  reads  this  number  of  equations  from  the  file. 

The  block  of  equations  in  the  core  are  then  partially  reduced,  resulting 
in  ones  at  the  locations  representing  the  diagonal  of  the  complete  matrix 
and  zeros  for  all  the  other  values  in  the  columns  with  ones.  This 
partial  reduction  of  the  equations  in  the  core  is  illustrated  schematically 
in  Figure  B.l  where  4  rows  of  a  small  16  x  16  matrix  are  shown. 

'  0  0  0  X  X  X  X  X  XX  XX  X  X  x" 
0100XXXXXXXXXXXX 
0010XXXXXXXXXXXX 
_0  0  0  1  X  X  X  X  X  X  X  X  X  X  X  X_ 

Figure  B.l  Illustration  of  How  a  Block  of  Equations  Represented  by 
a  Portion  of  a  16  x  16  Matrix  is  Partially  Reduced  in 
the  Core.  The  X's  Represent  Nonzero  Numbers. 


-108- 


26  June  1981 
AMY : cag 


With  the  operations  on  the  equations  in  the  core  completed,  the 

remaining  equations  or  rows  of  the  matrix  must  also  have  terms  eliminated 

by  the  equations  in  the  core.  Therefore,  the  next  procedure  is  to  read 

one  equations  at  a  time  from  the  file,  operate  on  it  with  each  equation 

in  the  core  and  write  the  remaining  terms  on  another  file.  After  all 

the  equations  have  been  operated  on  by  the  equations  in  the  block,  the 

first  step  of  the  reduction  is  complete.  After  the  first  step,  the  example 

16  x  16  A  matrix  would  appear  as  it  is  shown  schematically  in  Figure  B.2. 

As  with  the  other  reduction  technique,  the  operations  performed  on  the  A 

matrix  are  also  performed  on  the  onset  flow  matrices  which  will  represent 

the  solutions  when  the  reduction  is  complete. 

"lOOOXXXXXXXXXXX 
0100XXXXXXXXXXX 
0010XXXXXXXXXXX 
_0  0  0  1  XX  X  X  X  X  X  X  X  XX 

ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 
ooooxxxxxxxxxxx 

Figure  B.2  Schematic  of  a  16  x  16  Matrix  After  it  has  been 
Partially  Reduced  by  the  First  Four  Rows  of  the 
Matrix. 
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After  each  complete  step,  the  subroutine  reads  in  the  next  group  of 
equations  (or  rows  of  the  matrix)  and  repeats  the  reduction  process  until 
all  the  equations  have  been  in  the  reduction  block  in  the  core  and  the  matrix 
is  completely  reduced.  The  method  of  reducing  the  matrix  by  operating  on 
blocks  of  equations  which  are  in  the  core  of  the  computer  is  expected  to 
result  in  a  substantial  savings  in  computer  time  over  the  iterative  method, 
if  the  number  of  times  the  equations  must  be  read  from  file  is  less  for 
the  reduction  method.  The  number  of  times  the  equations  must  be  read 
from  a  file  is  determined  by  the  number  of  equations  which  will  fit  into 
the  allotted  space  in  the  core.  To  save  time,  the  subroutine  taken  advantage 
of  the  fact  that  it  need  not  operate  on  or  store  columns  of  the  matrix  which 
have  already  been  reduced  to  ones  and  zeroes.  Thus,  as  the  solution  proceeds 
and  the  number  of  terms  in  the  equations  is  decreasing,  a  larger  number  of 
equations  will  fit  into  the  allotted  area,  decreasing  the  required  number 
of  steps.  Figure  B.3  illustrates  the  second  step  of  the  reduction  for 
the  example  16  x  16  matrix  and  also  demonstrates  how  the  number  of  equations 
which  will  fit  into  the  core  increases  as  the  number  of  terms  decreases. 

To  better  realize  the  significance  of  increasing  the  number  of 
equations  in  the  core  as  the  number  of  terms  decreases,  a  more  realistic 
example  will  be  given.  Consider  the  case  where  there  are  200  equations 
each  with  200  terms  and  the  alloted  array  space  for  solving  the  equation 
is  50  x  50  or  2500.  Initially  200  spaces  are  required  for  reading  and 
writing  data  on  files,  thus  only  2300  spaces  are  available  for  storing  equations 
in  the  core.  With  200  terms  per  equations,  11  equations  can  be  stored  in  the 
core  and  11  columns  reduced.  If  the  solution  technique  did  not  take  advantage 
of  the  fact  that  after  each  step  the  number  of  terms  has  decreased,  it  would 
require  200  5  11  =  18.2  or  19  steps  to  reduce  the  entire  matrix.  However, 
with  the  number  of  terms  decreased  each  step,  more  equations  fit  into  the 
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Figure 
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B.3  Schematic  of  the  Example  16  x  16  Matrix  After 

the  Second  Step  of  the  Reduction.  The  Brackets 
Show  the  Block  of  the  Matrix  Which  was  in  the 
Core  for  the  Second  Step. 
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core  and  only  10  steps  are  required.  Table  B.l  depicts  the  steps  required 
for  this  example.  Ten  steps,  meaning  ten  times  the  matrix  must  be  read 
from  files  compared  to  approximately  20  iterations  required  for  the  iterative 
approach  would  result  in  a  substantial  savings  in  computer  time. 

One  potential  disadvantage  of  reducing  the  matrix  by  blocks  of  equations 
in  the  core  is  that  the  technique  does  not  utilize  the  concept  of  first 
arranging  the  matrix  for  a  strong  diagonal.  In  cases  run  to  date,  this 
has  not  been  found  to  cause  any  problems. 

Now  that  the  three  solution  techniques  and  their  core  storage  require¬ 
ments  have  been  discussed,  a  comparison  of  the  computation  time  for  each 
of  the  techniques  for  a  sample  problem  is  appropriate.  For  the  analytical 
airfoil  profile  in  a  linear  cascade  previously  discussed  in  the  section  on 
check  cases,  the  program  was  run  with  the  equations  solved  using  each  of 
the  three  techniques.  The  cascade  was  modeled  as  two  similar  cascades  to 
increase  the  number  of  points  and  to  check  the  multiple  cascade  solution 
procedure.  With  two  cascades  each  with  99  points  describing  the  airfoils, 
there  was  a  total  of  196  equations  each  with  196  terms.  Since  there  were 
two  cascades  with  circulation,  four  basic  solutions  were  required.  The 
reduction  technique  with  all  the  equations  in  the  core  utilized  an  array 
area  200  x  200.  The  reduction  of  the  matrix  with  only  a  block  of  equations 
in  the  core  was  performed  using  an  array  area  100  x  100.  Constructing  a 
table  similar  to  Table  B.l  would  show  that  the  matrix  was  reduced  in 
three  steps.  For  the  iterative  technique,  the  number  of  iterations 
required  by  each  of  the  basic  solutions  is  shown  in  Table  B.2  along  with  the 
computation  times  required  by  the  three  techniques.  Results  for  the 
iterative  technique  for  two  different  convergence  criteria  are  shown  in 


the  table. 
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TABLE  B.l 

Steps  in  Reducing  a  200  x  200  Matrix 
in  a  50  x  50  Array 


Step 

Terms  per 
Equation 

Array  Space 
Available  for 
Equations 

Equations  in 
the  Core  and 
Columns 
Reduced 

Total  Number 
of  Columns 
Reduced 

1 

200 

2300 

11 

11 

2 

189 

2311 

12 

23 

3 

177 

2323 

13 

36 

4 

164 

2336 

14 

50 

5 

150 

2350 

15 

65 

6 

135 

2365 

17 

82 

7 

118 

2382 

20 

102 

8 

98 

2402 

24 

126 

9 

74 

2426 

32 

158 

10 

42 

2458 

58  Maximum 

42  Needed 

200 

TABLE  B.2 

A  Comparison  of  the  Computation  Times  Required  by  the 
Three  Solution  Procedures  for  a  Sample  Problem 


Reduction 
with 
complete 
matrix  in 
the  core 

Reduction 
by  blocks 
of  equa¬ 
tions  in 
the  core 

Iterative  Procedure 

Convergence 
Criterion  =  .01 
Iterations  =  21, 
31,  21,  21 

Convergence 
Criterion  =  .001 
Iterations  =  26, 
41,  27.  26 

CPU  Time  (sec) 

232.55 

135.44 

102.58 

123.82 

I/O  Time  (sec) 

21.54 

31.03 

147.51 

189.50 

Total  Time  (sec) 

254.09 

166.47 

250.09 

313.32 
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It  is  seen  from  Table  B.2  that  the  procedure  which  reduces  the  matrix 
by  blocks  of  equations  in  the  core  required  the  least  amount  of  time  of 
the  three  techniques.  The  reduction  in  blocks  procedure  required 
significantly  less  CPU  time  than  the  procedure  with  the  complete  matrix 
in  the  core  for  two  reasons:  (1)  the  procedure  which  reduces  the  matrix 
by  blocks  does  not  even  consider  reduced  columns  where  as  the  other 
technique  merely  skips  zero  terms  when  encountered,  (2)  The  procedure  which 
reduces  the  complete  matrix  in  the  core  must  reorder  the  solution  when 
the  reduction  is  completed.  It  is  also  seen  from  Table  B.2  that  the 
reduction  in  blocks  technique  requires  considerably  less  computer  time 
than  the  iterative  technique  because  of  its  savings  in  I/O  time. 

I/O  (i.e.  input/output)  time  reflects  the  time  spent  reading  disc 
files.  The  iterative  technique  must  read  the  complete  set  of  equations 
from  a  file  for  each  iteration;  and  as  shown  in  the  table,  for  a 
convergence  criteria  of  .001  the  second  basic  solution  required  41 
iterations. 

To  conclude  this  section  of  Appendix  B,  Table  B.3  summarizes  the 
array  storage  requirements  and  the  major  advantages  and  disadvantages 
of  the  three  solution  procedures.  With  a  general  understanding  of  the 
three  procedures  and  there  storage  requirements  from  this  section,  the 
following  section  explains  specifically  the  array  dimensions  and  how 
they  and  the  solution  procedure  can  be  changed. 
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TABLE  B.3 

A  Comparison  of  the  Three  Solution  Procedures 


Solution 

Procedure 

Array  Storage 
Spaces  Required 

Major 

Advantages 

Major 

Disadvantages 

Reduction  with  the 
complete  matrix  in 
the  core 

N2 

Utilizes  strong 
diagonal  concept 

1.  Limited  to  small 
problems  because  of 
storage  requirements 

2.  Longer  computa¬ 
tion  time 

Reduction  by  blocks 
of  equations  in  the 
core 

At  least  2N,  larger 
storage  areas  will 
reduce  the  compu¬ 
tation  time 

1.  Shorter  computa¬ 
tion  time 

2.  Can  utilize  what¬ 
ever  storage  space 
is  available 

Does  not  utilize 
strong  diagonal 
concept 

Iterative 

Procedure 

N 

Small  array 
storage  require¬ 
ments 

Longer  computation 
time 

N  =  number  of  equations  (and  number  of  terms  per  equation) 
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B. 2  Changing  Array  Dimensions  and  Solution  Procedures 

It  may  sometimes  be  necessary  to  change  the  array  dimensions  to  match 
the  program  to  a  particular  computer  system  or  a  particular  problem. 
Changing  solution  procedures  primarily  involves  making  certain  the  array 
dimensions  are  consistent  with  the  procedure,  thus  both  changing  array 
dimensions  and  solution  procedures  will  be  discussed  in  the  context  of 
explaining  the  arrays  in  the  program.  Several  segments  of  the  computer 
program  are  presented  in  figures  to  aid  in  the  discussion.  All  the  major 
arrays  used  in  the  program  are  in  common  statements,  utilizing  the  same 
storage  space  for  all  the  subroutines  and  transferring  the  needed  values 
from  subroutine  to  subroutine. 

Figure  B.4  is  a  listing  of  the  main  program  of  the  rewritten  version 
of  the  Douglas  Neumann  program  which  reveals  all  the  major  arrays  used. 

The  COMMON  statements  BLK2  and  BLK3  contain  tne  arrays  which  describe 
the  geometry  of  the  airfoils  in  the  cascades.  The  variable  arrays  in 
these  statements  are  defined  below: 

X,  Y  -  the  coordinates  input  to  the  program  describing  the  shape 
of  the  airfoils 

XMP,  YMP  -  the  coordinates  of  the  midpoints  of  the  surface  elements 
defined  by  X,  Y 

Q  -  complex  body  coordinates  (X  +  iY) 

SINA,  COSA  -  sine  and  cosine  of  the  angle  of  the  surface  elements 

Since  the  arrays  in  the  COMMON  statements  BLK2  and  BLK3  contain 
airfoil  coordinate  data,  it  is  obvious  that  they  must  be  dimensioned 
as  large  as  the  number  of  points  used  to  describe  the  shape  of  the  airfoils 
in  the  cascades.  The  current  version  of  the  program  was  intended  to 
handle  up  to  10  different  cascades  each  of  which  can  usually  be  adequately 
described  with  100  points  or  less.  These  numbers  (10  x  100)  resulted  in 
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C0MM0N/BLK1 /HEDRI20) .THETA 

_  COMMON/BLK2/  Q(  10.00). SINAU000)  .CQSAC  1000) _  ._ _ . _ 

COMMON/BLK3/  X(1000).Y(1000 ) ,XMP( 1000 ) .  YMP( 1000 ) 

C  C0MM0N/3LK2/  O (  3 0 0  )  , S I NA ( 3 0 0 ) . COS A < 30 0 ) 

C _ COMMON/ BLK 3/  X(30 0  )  •.Y(_30 0  >_.*  XMP  (  30  0.)  t  YMP(  30_0  ) _ _ _ 

C0MM0N/8LK4/  NP( l 0) *NLE( 1 0 > . SUMDS { 1 0 > . CHOR D< 10) .XP (10),YP{ 10) • 
*RPADDX(  10) . CC  ANG (  1O).AD0Y( 10).NLF{  10) 


5  _ C.OMMQN^RLK5/  FLG0  2.tFLG03  FL  G04  ,  FLG05.  Fl_ G06  .  FLGO  7.  FLG08  .  F.LG09  . 

♦FLG10 .FLG1 1 .FLG12 

6  COMMON/8LK6/  NB . NT . I M . NS DL . RPI . R2 P 1 » S® . CL . ALPHA . F ALPHA , O ALPHA 

Z _ C.QM  MON/  f)  L  <  7  /.  XQBC500)  ,.YOB(.500  )  .  NO  BP...  _ _ _ 

S  C0MM0N/BLK8/  SIG< 1000.12) 

C  COMMON/QLKfl/  SIG{300.12) 

9 _ _  C0MM0N/BLK9/  A(  40000)  .  . . . . . . . . 

C  ARRAY  A  IN  COMMON  BLOCK  9  IS  USED  BY  SEVERAL  DIFFERENT 

C  SUBROUTINES  TO  STORE  TEMPORARY  VARIABLES 

1.3 _ COMPLEX..! M.Q  _ ; _ 

11  INTEGER  FLGO? .FLG03. FLGO  4  «  FLGO 5. FLG06 . FLGO 7 . FLG08. FLG09* 

♦FLG10.FLG11 .FLG12.SU8KS 


C 

c 

_ _ c _ 

1  ?  10 
1  3 

_ _ c 

c 

c 

1  4 
15 

_ c _ 


PARTI — READS  IN  OATA,  SCALES,  ROTATES  AND  TRANSLATES  THE  AIRFOILS. 
IT  ALSO  TRANSFORMS  RADIAL  FLOW  CASCADES 


CALL  PARTI 

IF ( FLG08.NE  .0  )  GO  TO  30 


PART2 — SETS  UP  THE  MATRIX  WHICH  IS  SOLVED  FOR  THE  SOURCE 
DISTRIBUTIONS 


CALL  PART  2 

IFINT .LE.200)  GO  TO  20 


C  3ART3 — SOLVES  THE  MATRIX  BY  REDUCING  THE  EQUATIONS  IN  BLOCKS 

C  IF  THE  NUMBER  OF  ELEMENTS  IS  GREATER  THAN  200 

16  CALL  PART3 

17  GO  TO  30 

C  PAR T 4  —  SOLVES  THE  MATRIX  BY  REDUCTION  IF  THE  NUMBER  OF  ELEMENTS 

C  IS  LESS  THAN  OR  EQUAL  TO  200 

13  20  CALL  PART  4 

C 

C  PAPT5  —  SOLVES  THE  COMBINATION  EQUATIONS*  CALCULATES  THE  VELOCITIES. 

C  AND  PRESSURES  AND  WRITES  THE  OUTPUT 

C 

19  . . .  30.  CALL  PART 5. _ _ _ _ _ _  _ _ 

20  GO  TO  10 

21  END 


Figure  B.4  A  Listing  of  the  Main  Program 
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the  dimension  of  100  seen  in  Figure  B.4.  Any  number  of  points  can  be  used 
to  describe  the  airfoils  of  cascades,  provided  the  total  does  not  exceed 
the  dimensions  in  COMMON  statements  BLK2  and  BLK3.  For  many  problems, 
however,  a  large  number  of  points  is  not  required  and  the  dimensions  in 
COMMON  statements  BLK2  and  BLK3  can  be  reduced.  Shown  as  comment  statements, 
the  COMMON  statements  with  the  arrays  dimensioned  300  are  an  example  of 
the  dimensions  used  for  several  sample  runs.  This  example  is  carried 
throughout  the  other  COMMON  statements  and  subroutines  to  illustrate  the 
changes  which  must  be  made  to  keep  all  the  array  dimensions  consistent. 

COMMON  statement  BLK4  contains  various  cascade  parameters.  Since 
the  program  was  intended  to  handle  up  to  10  cascades,  the  variables  in 
COMMON /BLK.4 /  are  all  dimensioned  10.  These  dimensions  could  be  reduced 
for  problems  with  fewer  than  10  cascades.  However,  reducing  these 
dimensions  would  yield  an  insignificant  savings  in  storage  space  and  it 
is  not  recommended  they  be  reduced.  If  it  every  is  desirable  to  analyze 
more  than  10  cascades,  the  dimension  of  the  variables  in  C0MM0N/BLK4/  would 
have  to  be  increased. 

The  next  major  COMMON  statement  is  C0MM0N/BLK7/,  which  contains  the 
coordinates  (XOB,  YOB)  of  the  off-body  calculation  points.  The  off-body 
point  coordinates  are  stored  in  separate  arrays  from  the  body  points  in  the 
rewritten  version  of  the  Douglas  Neumann  program  so  that  the  number  of 
off-body  calculation  points  can  be  changed  easily  without  affecting  the 
number  of  body  points. 

The  statement  C0MM0N/BLK8/  SIG(1000,12)  dimensions  the  SIG  array  which 
stores  the  distribution  of  sources  on  the  surfaces  of  the  airfoils.  Each 
surface  element  has  a  constant  value  of  source  strength  for  each  of  the 
basic  solutions  required.  The  first  dimension  of  SIG  must  therefore  be 
greater  than  or  equal  to  the  number  of  surface  elements,  and  the  second 
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dimension  must  be  greater  than  or  equal  to  the  number  of  basic  solutions.  If 
there  are  10  cascades  with  circulation  to  be  analyzed,  there  are  10  +  2  =  12 
basic  solutions  required.  SIG  is  therefore  dimensioned  (1000,  12)  in  Figure  B.4 
which  is  consistent  with  the  previous  dimensions  in  BLK2,  BLK3  and  BLK4.  In  the 
comment  common  statement,  SIG  is  dimension  (300,  12),  to  be  consistent  with 
comment  C0MM0N/BLK3/  and  C0MM0N/BLK4/  shown  in  the  figure.  The  300  dimension 
in  BLK3  and  the  12  dimension  is  unchanged  since  the  dimensions  in  BLK4  are 
unchanged.  In  short,  the  SIG  dimensions  can  be  changed,  but  they  should  be 
consistent  with  the  dimensions  in  the  other  arrays. 

The  final  COMMON  storage  area  is  COMMON/ BLK9 / .  The  A  array  is 
COMMON/BLK9/  is  used  to  solve  the  set  of  simultaneous  equations.  In 
SUBROUTINE  PARTS  the  same  common  space  is  divided  among  three  arrays 
which  are  used  to  calculate  the  velocities  on  the  airfoils  and  reorder  the 
data.  Since  C0MM0N/BLK9/  is  dimensioned  several  different  ways  in  the 
various  subroutines,  care  must  be  taken  to  make  all  changes  consistent 
throughout  the  program. 

The  program  normally  contains  two  methods  of  solving  the  set  of 
simultaneous  equations  and  selects  one  method  according  to  the  size  of 
the  problem  being  run.  The  equations  are  either  solved  by  reduction  of 
the  complete  matrix  in  an  array  in  SUBROUTINE  PART4  and  SUBROUTINE  MIS1,  or 
the  equations  are  solved  by  reduction  of  blocks  of  equations  or  by  iterating 
in  one  of  the  two  versions  of  SUBROUTINE  PART3 .  The  program  user  can  inter¬ 
change  the  two  versions  of  SUBROUTINE  PART3  as  he  chooses.  The  dimension 
of  the  A  array  determines  the  maximum  number  of  equations  which  can  be 
solved  using  the  straight  matrix  reduction.  In  Figure  B.4,  the  A  array 
is  dimensioned  40000  which  corresponds  to  a  square  matrix  200  x  200.  Thus, 
with  this  dimension,  up  to  200  equations  can  be  solved  by  reducing  the  complete 
matrix  in  an  array.  The  statement  IF(NT.  LE.  200)  go  to  20  shown  in  Figure  B.4 
directs  the  program  to  SUBROUTINE  PART4  if  the  number  of  elements  (equations) 
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is  less  than  or  equal  to  200.  The  number  in  the  IF  statement  may  be 
changed  to  be  less  than  the  square  root  of  the  A  dimension,  but  it  may 
not  exceed  the  square  root  of  the  A  dimension. 

Figure  B.5  shows  the  beginning  segments  of  subroutines  PART4  and 
MIS1  and  reveals  the  COMMON  statements  required  by  these  subroutines. 

It  should  be  noted  that  in  SUBROUTINE  PART4  the  A  array  is  two-dimensional 
and  has  equal  dimensions  since  it  represents  a  square  matrix.  Because 
of  its  form  in  SUBROUTINE  PART4,  it  is  clear  that  the  dimension  of  the 
A  array  in  other  subroutines  must  be  such  that  it  can  be  represented  as 
a  square  two-dimensional  array.  It  is  also  shown  in  Figure  B.5  that  in 
SUBROUTINE  MIS1  both  the  A  array  in  C0MM0N/BLK9/  and  the  B  array  in 
C0MM0N/BLK8/  are  one-dimensional  while  in  the  calling  subroutine  (SUBROUTINE 
PART4)  both  arrays  are  two-dimensional.  The  parameters  ND  and  ND2  in 
SUBROUTINE  MIS1  provide  MIS1  with  the  first  dimension  of  A  and  B  in  their 
two-dimensional  form  and  enable  MIS1  to  recover  values  from  the  one-dimensional 
arrays  which  were  stored  as  values  in  two-dimensional  arrays.  If  the 
dimensions  of  A  or  B  are  changed,  the  parameters  ND  and  ND2  must  also 
be  changed. 

C0MM0N/BLK9/  also  appears  in  the  two  versions  of  SUBROUTINE  PART3 . 

In  the  version  which  solves  the  equations  by  reducing  blocks  of  equations, 
the  A  arrays  scores  the  equations  being  reduced.  The  beginning  segment 
of  this  version  of  SUBROUTINE  PART3  is  shown  in  Figure  B.6.  It  is  seen 
from  the  figure  that  the  A  array  is  one-dimensional.  The  parameter 
NA  shown  in  Figure  B.6  must  be  given  the  value  of  the  A  array  dimension 
in  order  that  the  subroutine  can  calculate  the  number  of  equations  which 
will  fit  into  the  array.  Thus,  if  the  dimension  of  A  is  changed,  NA 
must  also  be  changed.  The  same  COMMON  statements  are  also  used  by  the 
version  of  SUBROUTINE  PART3  which  iterates  for  a  solution.  No  parameters 
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596  SUBROUTINE  PART  4 
C 

_ C _ IP  THE  NUMBER  OF  ELEMENTS  IS  LESS  THAN  OR  EQUAL  _TO.  200  *  THI  S 

C  SUBROUTINE  SETS  UP  THE  MATRICIES  IN  ARRAYS  WHICH  ARE  THEN  SOLVED 

C  BY  REDUCTION  IN  SUBROUTINE  MISi 

.  .C  ....  . .  ...... 

597  C  OMMON/RLK  2 /  Q (  1 0 0 0  )  . S IN A {  1 0 0 0  ) . COSA (  l 0 00 ) 

C  C  OMMON/BLK  2/  Q {  3 0 0 ) , S I NA <  3 0 0 ) . CO S A ( 30 0 ) 

595 _  _  COMMQN/OLK6/  NB,  NT  *  I M.  NS  CL  .  RPI  ,R2PI  *  S°.CL  •  ALPHA  ,  F  ALPHA  ,  O  ALPHA 

599  COMMON /BLK3/  SIG<I000.12) 

C  C0MM0N/BLK8/  SIGI300.12) 

600  C0MMON/BLK9/  A(200.200)  _  . . . 

601  COMPLEX  IM.Q 


623  SUBROUTINE  M I S 1 ( N . NSOL . NE RR ) 

C  THIS  SUBROUTINE  REDUCES  THE  A  MATRIX  FOR  EACH  VECTOR  IN  B. 

_ C _ .''THE..  FINAL  SOLUTIONS  ARE  IN  8« _ _ _ 

624  A0MMON/BLK8/  8< 12000) 

C  dfJMMON/BLKS/  3(  36001 

625  _ _ C ^IMMON/ BLK9/.  A {  400  00  ) _ _ _  _ _ _ _ _ _ _ 

C  CHANGE  N02  WHEN  CHANGING  THE  DIMENSION  OF  B 

C  NCi2  =  (B  D  I  MENS  I  ON  )  /  l  2 

_ C _ NO  s _ S ORTIA  DIMENSION  ) _ _ _  _ _ _ 

C 

626  NERR=1 

627  ....  ND  =  200  . . . . . . .  . . 

625  N02=1000 

629  M=NSOL 


Figure  B.5  Parts  of  Subroutines  PART4  and  MISI  Showing  the 
Array  Dimensions 
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456  SUBROUTINE  PART3 
C 

. . C _  IF  . .THE.  NUMBER  .OF  SQUAT  I  ONS  I  S.  GREATER.  THAN  .  200  i _ 

C  THIS  SUBROUTINE  SOLVES  THE  EQUATIONS  BY  REDUCTION 

C  HAVING  ONLY  A  PARTIAL  BLOCK  OF  EQUATIONS  IN  THE  CORE 

_ C. . AT  ONE  T IMS.  THE  NUMBER  OF  EQUATIONS  IN  THE_ BLOCK  INCREASES _ 

C  AS  THE  NUMBER  OF  TERMS  IS  REOUCED,  THUS  DECREASING  THE 

C  STEPS  RFOUIREO. 

457  COMMON/BLK 2 /  Q <  1 00 0 ) • S  INA (  1 0 0 0 ) * COSA C  1 0 00 ) 

C  C0MM0N/3LK2/  Q ( 3 0 0  )  .  S I NA ( 3 0 0 > • COS A ( 3 0 0 > 

45  3  C0MM0N/BLK6/  _NB.NT  ,  IM.NSOL.RPI  . R2P I . SP *  CL . ALPHA , F ALPHA, D ALPHA 

450  C0MM0N/3LK9/  SIGI1000.12) 

C  C0MM0N/8LK8/  SIG(300.12> 

_4.60. _ C.aMMQN./BLK9/_..A  i  4.0.0  0  0  )_ _ _ : _ 

461  COMPLEX  IM.Q 

462  REWIND  3 

463  _  REWIND..  B . . . . _ _ _  _  _ _ _ _ _ 

464  WRITEI6.601 ) 

465  601  =0RMAT ( • 0 •  ,  10 X . • THE  EOUATIONS  ARE  REDUCED  IN  BLOCKS  ', 

_ _ * _ _ _ _ •  IN.  SUQROUT.I  NE..P  AR  T  3.*. )... _ 

466  DO  30  1=1, NT 

467  S I G ( I  *  1  ) =S I NA {  I  ) 

469  _  _ _ SIG(  l  ,2)  =-COS  AC  I  )  _  _ 

460  I = ( NSOL.LT. 3)  GO  TO  30 

470  DO  10  J=3,NS0L 

47.1 _ REA01.3I  S.JLGLIjlJJL _  _  _ _ _ 

472  10  CONTINUE 

473  DO  20  J=3.NSQL 

474  .  . . SIG<  I  ,  J)=— SIG(  I  ,  J)  . .  .... 

475  20  CONTINUE 

476  30  CONTINUE 

C  NEQB  =  NUMBER  OF  EQUATIONS  IN  THE  BLOCK,  NEQB  CHANGES 

C  AS  THE  NUMBER  OF  TERMS  IS  REDUCED. - 

. . _.C _ NJNR.F  NUMBER  .OF  TERMS.  I  N  .THE  .EQUATIONS  WHICH  HAVE  NOT _ _ _ 

C  BEEN  REDUCED  TO  0  OR  I 

C  NTR  =  NUMBER  OF  TERMS  PER  EQUATION  REDUCED  TO  0  OR  1 

..._ _ C _ NA  =  DIMENSION  OF  THE  A  ARRAY.  _ _ _ 

C 

C  THE  LAST  SEGMENT  OF  THE  A  ARRAY  IS  USED  FOR  TEMPORARY 

_C  STORAGE  REQUIRED  FOR. READING  AND  WRITING  ON  FILE. . . 

C  IFILE  =  CODE  FOR  DETERMINING  WHICH  FILE  IS  READ  AND  WHICH 

C  IS  WRITTEN  ON 

477  NA  =  40  00  0 

479  I F I LE  =  8 

479  _  NTR=0  .  _  _  ,  _ _ _  _ 

430  40  NTNR=  NT -NTR 

491  NE  08=  <  N A— NT  NR ) /NT  NR 

492  IFINEQB.GT.NT)  NEQB=NT.  _  _ 

Figure  B.6  Beginning  at  SUBROUTINE  PART3  Showing  the  COMMON  Statements 
and  Important  Parameters 
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need  to  be  changed  in  the  iterative  version  when  array  dimensions  are 
changed. 

When  the  set  of  simultaneous  equations  have  been  solved  for  the  source  distri¬ 
butions  on  the  airfoils,  the  storage  space  in  C0MM0N/BLK9/  is  free  and  can  be  used 
for  other  purposes.  This  area  is  therefore  used  in  SUBROUTINE  PART5  for  calcu¬ 
lating  velocities  and  rearranging  the  data.  The  beginning  segment  of  SUBROUTINE 
PART5  shown  in  Figure  B.7  reveals  that  the  storage  space  in  C0MM0N/BLK9/  is  divided 
among  three  arrays.  The  space  can  be  divided  in  any  manner  among  the  array,  provided 
each  array  has  sufficient  space  to  perform  its  task.  The  B  array  is  used  to 
store  the  coefficients  for  the  influence  of  the  sources  of  all  the  elements 
on  a  particular  element.  The  B  matrix  is  very  similar  to  the  A  matrix  of 
Equation  B.l  except  the  velocity  component  calculated  from  the  B  matrix  is 
tangent  to  the  body  surface  rather  than  normal  to  the  surface.  Only  one 
row  of  the  complete  B  matrix  needs  to  be  stored  in  the  B  array  at  one  time, 
thus  its  dimension  can  be  as  small  as  the  total  number  of  surface  elements. 

As  the  velocities  at  the  midpoints  of  the  elements  are  calculated,  the  values 
are  stored  in  the  array  V.  The  velocities  are  calculated  for  one  cascade 
and  then  the  results  are  printed.  Therefore,  the  V  array  need  only  store 
the  velocities  for  one  cascade  at  a  time  and  V  must  be  dimensioned  at  least  as 
as  large  as  the  largest  number  of  elements  representing  any  airfoil  in  an 
analysis.  The  array  TSTOR  is  used  to  temporarily  store  the  velocities  and 
element  midpoint  coordinates  for  the  cases  where  the  data  must  be  reordered 
at  the  end  of  the  computations.  TSTOR  must  also  be  dimensioned  at  least 
large  as  the  largest  number  of  elements  representing  any  airfoil.  For 
example,  in  a  multiple  cascade  analysis  with  three  cascades  represented 
by  80,  90,  and  100  elements,  V  and  TSTOR  must  both  be  dimensioned  at  least 
100.  For  this  example  B  must  be  dimensioned  at  least  as  large  as  the  total 
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637  SUBROUT INS  PARTS 

C 

....  C _  THIS.  SUSKO.UTI  NE...SOLVES  THE..  COMB  I  NAT  I  ON  EQUATIONS  ..AND  _ .'. _ 

C  CALCULATES  the  velocities  and  pressures. 

c 

683 _ _ _ COMMON/ J3LK  I  /HEDR (  2 0  )  »  THET  A  _ _ _ _  _  .  _ _ _ 

639  COMMON/BLK2/  0 ! 1 0 0 0 )  , S I N A !  I  0 0 0  > , COS A ! 1 0 00 > 

690  C0MM0N/9LK3/  X (  1 0 0 0 ) . Y < 1 0 00  )  . X MP (  1 00 0 )  » YMP (  1  000 ) 

_  . jC _ COMMON/BL<2./  Q  (  30  0  )  .  S  I NA  (  300  )  .  CQS  A  (  300  ) _ _ 

C  COMMON/BLK3/  X < 3 0 0 ) . Y{ 30 0 ) , XMP ( 3 0 0 ) . YMP( 30 0 ) 

691  C0MM0N/8LK4/  NP< 10 ) • NLE( 101.SUMDS! 10 » .CHORD! 1 0 ) . XP { 1 0 ) , VP ( 1 0 ) 

*R»ADDX 1  10 )  .CCANG! I  0  )  , AODY ! 1 0) . NLF! 10)  .  .  . . 

6  92  COMMON/ 3LK5/  EL G02 • FLG03 . FLG04  * FLG 05 * FLG06 • FLGO 7 . F LGO 8 . FLG 09 . 

♦FLG10 .FLG1 1 .FLG 1 2 

69  3  ...  ....  _  COMMON/BuKS/  NB.,.NT.  I  M»NS.aL.».R=,I_».R2PI.».S?..CLjLALPHA,F  ALPHA,  D  ALPHA 

694  COMMON/9LK7/  XOS! 500). YOB ( 500 ) , N03P 

695  COMMON/BLK8/  SIG(1000.12) 

...  _ C  C0MM0N/3LK9/.  SIG(  300, 12)  ....  _  _ _ 

695  C0MM0N/BLK9/  3 <  l 60 0 0  )  . V (  I  200 0 ) . TST OR (  1 20 0 0 ) 

697  COMPLEX  I M , Q, CSI NH, ZOB. TF , W 

693  _  ._I  NTE3FR  FLG02  ,  FLGO  3  .  FLG04  ,  FL  GO  5 ,  FLGO  6  ,  FL  G07  ,.FLG  08  ,  FLG  O.9., 

4FLG10.FLG1 1 , F  LG  I  2 , SUBKS 

699  DIMENSION  O VB AS <  I  0 , 1 2 ) * D V( 1 1 , 1 1 ) • DVS ( I  1 . I  1 ) * GAM!  1 1  ) 

.700  DIMENSION  COMBS!  12)  _  . 


Figure  B.7  Beginning  of  SUBROUTINE  PART5  Showing  the 
Array  Dimensions 
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number  of  elements  or  80  +  90  +  100  =  270.  Under  usual  circumstances 
C0MM0N/BLK9/  will  contain  much  more  space  than  actually  required  by 
B,  V  and  TSTOR. 

The  major  arrays  in  the  program  have  now  been  discussed  and,  if  it  is 
necessary,  the  program  user  should  be  able  to  change  the  dimensions  of 
the  arrays.  The  user  can  also  change  the  criterion  which  directs  the 
program  to  one  of  the  two  internal  solution  procedures  or  interchange 
the  two  versions  of  SUBROUTINE  PART3,  since  the  primary  concern  when 
making  such  changes  is  providing  adequate  and  consistent  array  storage. 
When  changing  the  dimensions  of  arrays  which  are  used  for  solving  the 
set  of  simultaneous  equations,  the  user  must  remember  to  change  the 
appropriate  parameters  in  the  subroutines  as  discussed  in  this  appendix. 
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APPENDIX  C 

Coordinate  Data  Input  Order  and  Reordering 

The  coordinate  data  describing  the  shape  of  the  bodies  in  a  cascade 
are  obviously  important  input  to  the  Douglas  Neumann  program.  Of  great 
importance,  although  for  reasons  not  as  obvious,  are  the  starting  point 
and  order  of  the  data.  The  order  of  the  data  is  important  because  the 
sine  and  cosine  of  the  angle  of  the  surface  elements  are  calculated  from 
the  coordinate  data.  If  the  coordinate  data  are  not  in  the  correct  order, 
the  sign  of  the  angle  according  to  the  chosen  convention  will  be  incorrect, 
and  the  boundary  condition  will  not  be  properly  satisfied  on  the  body 
surface.  For  linear  cascades  and  circular  cascades  with  the  flow  radially 
inward,  the  starting  point  is  also  important  because  it  tells  the  location 
of  the  trailing  edge.  The  location  of  the  trailing  edge  must  be  known  for  the 
program  to  satisfy  the  Kutta  condition. 

Concerning  the  starting  point  and  order  of  the  body  coordinate  data, 
there  are  three  types  of  cascades  which  must  be  given  separate  considera¬ 
tion.  The  three  types  are  the  linear  cascade,  the  circular  cascade  with 
the  radial  flow  inward,  and  the  circular  cascade  with  the  radial  flow 
outward.  To  avoid  confusion,  it  is  desirably  to  keep  the  starting  point 
and  order  of  the  data  consistent  for  all  three  types  of  cascades.  The 
order  of  the  data  was  chosen  to  match  the  original  Douglas  Neumann  program 
for  a  linear  cascade.  For  the  circular  cascades,  the  program  must  reorder 
the  data  so  that  after  the  transformation,  the  starting  point  and  order 
correspond  to  the  convention  for  the  linear  cascade.  In  this  Appendix, 
the  coordinate  data  input  for  each  type  of  cascade  will  be  described  and 
the  reordering  of  the  data  performed  by  the  program  will  be  discussed. 


Jk 
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For  all  three  types  of  cascades,  the  coordinate  data  describing  the 
shape  of  a  body  must  be  input  to  the  program  in  terms  of  a  local  coordinate 
system  x'-y'.  An  example  of  an  airfoil  in  its  local  x’-y'  coordinate 
system  is  shown  in  Figure  Cl.  This  airfoil  will  be  used  for  each  type 
of  cascade  to  illustrate  how  the  local  coordinate  system  for  the  airfoil  is 
related  to  the  cascade  coordinate  system.  As  illustrated  in  Figure  Cl, 
the  coordinate  data  must  start  at  the  trailing  edge  and  progress  clockwise 
around  the  airfoil.  (Note:  For  circular  cascades,  the  orientation  of 
the  x'-y'  coordinate  system  with  respect  to  the  cascade  coordinate  system 
is  the  same  regardless  of  whether  the  flow  is  radially  inward  or  outward. 

Thus,  the  true  leading  edge  and  trailing  edge  may  be  reversed  from  the 
orientation  shown  in  Figure  Cl.  Subsequent  figures  will  clarify  this 
point. ) 

Figure  C2  shows  how  the  airfoil  and  its  local  coordinate  system  relate 
to  the  cascade  coordinate  system  x-y  for  the  linear  type  cascade.  Although 
the  local  coordinate  system  of  the  airfoil  can  be  rotated  and  translated 
within  the  global  system,  the  cascade  must  parallel  the  y  axis  and  the  leading 
edge  of  the  airfoil  must  be  toward  the  -x  direction.  The  component  of  the 
flow  perpendicular  to  the  cascade  is  from  left  to  right  as  shown  in  Figure 
C2.  Since  the  cascade  solution  technique  is  the  same  for  all  three  types 
of  cascades,  the  cascade  shown  in  Figure  C2  also  represents  the  orientation 
of  the  cascade,  the  direction  of  flow,  and  the  order  and  starting  point 
of  the  data  for  the  radial  flow  cascades  after  they  are  transformed. 

The  orientation  of  the  local  coordinate  system  for  an  airfoil  in  a 
circular  cascade  with  the  radial  flow  direction  inward  is  shown  in  Figure  C3-a. 
Again  the  airfoils  can  be  rotated  and  located  at  any  R-0  position,  but  the 
origin  of  the  x'-y’  system  must  be  facing  R  ■  After  the  cascade  is 
transformed  to  a  linear  cascade,  the  orientation  of  the  airfoils,  the 
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L.E. 


Figure  C. 


Y' 


DIRECTION  FOR  COORDINATE  INPUT, 
STARTING  AT  THE  T.E. 


.  An  Airfoil  in  its  Local  Coordinate  System  x'-y',  With 
the  Order  for  the  Data  Input  Indicated  by  the  Arrows 
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Figure  C.2.  A  Linear  Cascade  Showing  the  Orientation  of  Airfoil 

Coordinate  System  With  Respect  to  the  Cascade  Coordinate 
System 
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RADIAL  FLOW 
DIRECTION 
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Figure  C.3.  a.  Schematic  of  a  Circular  Cascade  Showing  The  Orientation  of 

the  Airfoil  Coordinate  System  in  the  Global  System,  b.  A  Schematic 
Showing  the  Cascade  Orientation  and  Coordinate  Input  Direction 
After  the  Transformation,  c.  Final  Cascade  Configuration  and 
Order  of  Coordinate  Data 
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order  of  the  coordinate  data  and  the  flow  direction  are  shown  schematically 

in  Figure  C3-b.  Figure  C3-b  reveals  that  the  V  direction  does  not 

x 

correspond  to  the  -x  to  +x  direction  required  by  the  cascade  solution 

procedure.  For  this  reason,  an  additional  transformation  step  is  performed 

to  change  the  sign  of  the  x  coordinates  of  the  transformed  cascade, 

resulting  in  the  cascade  shown  in  Figure  C3-c.  After  this  step.  Figure 

C3-c  shows  that  the  cascade  orientation  and  V  flow  direction  are  correct, 

x 

but  it  is  seen  that  the  coordinate  data  as  read  in  do  not  correspond 
to  the  order  choosen  for  the  cascade  solution  procedure.  Rather  than 
having  different  input  orders  for  different  types  of  cascades,  the 
computer  program  reverses  the  order  of  the  coordinate  data  and  the 
cascade  solution  proceeds.  Following  the  flow  field  solution,  the  order 
of  the  data  is  again  reversed  so  the  output  order  matches  the  original 
input. 

Figure  C4-a  shows  the  orientation  of  the  x'-y*  system  for  an  airfoil 
in  the  third  type  of  cascade,  a  circular  cascade  with  the  radial  flow 
outward.  The  x'-y'  orientation,  order  of  the  coordinate  data,  and  the 
starting  point  are  the  same  as  for  a  circular  cascade  with  the  radial 
flow  inward.  Identical  coordinate  input  for  the  two  types  of  circular 
cascades  was  choosen  for  the  computer  program  so  that  wicket  gates  and 
stay  vanes  of  hydraulic  pump-turbines  could  be  analyzed  in  both  modes  with 
the  same  coordinate  data.  A  code  in  the  input  tells  the  program  the  direction 
of  the  flow,  and  for  the  circular  cascade  with  the  radial  flow  outward 
the  location  of  the  trailing  edge  is  specified  since  in  this  case  the  data 
do  not  start  there. 

When  the  cascade  in  Figure  C4-a  is  transformed,  the  resulting  linear 
cascade  is  shown  schematically  in  Figure  C4-b.  It  is  seen  in  this  figure 
that  the  orientation  of  the  airfoils  and  the  direction  of  V  are  correct 

V 


-131- 


26  June  1981 
AMYicag 


for  the  cascade  analysis.  However,  the  data  start  at  the  point  which 
would  be  the  trailing  edge  if  the  flow  was  radially  inward,  and  the 
starting  point  must  be  moved  to  the  true  trailing  edge  for  the  radial 
flow  outward  condition.  Again  the  computer  program  reorders  the  data 
so  that  the  input  can  be  consistent  for  each  type  of  cascade.  For  the 
circular  cascade  with  the  radial  flow  outward,  the  index  of  the  coordinate 
pair  which  give  the  location  of  the  trailing  edge  is  input  to  the  program, 
and  this  point  is  made  the  starting  and  finishing  point.  As  in  the 
case  when  the  radial  flow  is  inward,  the  output  is  put  back  into  the 
original  order  after  the  flow  field  calculations  are  complete. 


DIRECTION  OF 
COORDINATE  INPUT 


Figure  C.4.  a.  Schematic  of  a  Circular  Cascade  with  the  Radial  Flow  Direction 
Outward.  The  x'-y*  Orientation  is  the  Same  as  in  Figure  C.3a. 
b.  Final  Configuration  of  the  Cascade  After  the  Transformation 
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